Chapter 6 Diversity analysis

6.1 Alpha diversity

# Calculate Hill numbers
richness <- genome_counts_filt %>%
  column_to_rownames(var = "genome") %>%
  dplyr::select(where(~ !all(. == 0))) %>%
  hilldiv(., q = 0) %>%
  t() %>%
  as.data.frame() %>%
  dplyr::rename(richness = 1) %>%
  rownames_to_column(var = "sample")

neutral <- genome_counts_filt %>%
  column_to_rownames(var = "genome") %>%
  dplyr::select(where(~ !all(. == 0))) %>%
  hilldiv(., q = 1) %>%
  t() %>%
  as.data.frame() %>%
  dplyr::rename(neutral = 1) %>%
  rownames_to_column(var = "sample")

phylogenetic <- genome_counts_filt %>%
  column_to_rownames(var = "genome") %>%
  dplyr::select(where(~ !all(. == 0))) %>%
  hilldiv(., q = 1, tree = genome_tree) %>%
  t() %>%
  as.data.frame() %>%
  dplyr::rename(phylogenetic = 1) %>%
  rownames_to_column(var = "sample")

# Aggregate basal GIFT into elements
dist <- genome_gifts %>%
  to.elements(., GIFT_db2) %>%
  traits2dist(., method = "gower")

functional <- genome_counts_filt %>%
  column_to_rownames(var = "genome") %>%
  dplyr::select(where(~ !all(. == 0))) %>%
  hilldiv(., q = 1, dist = dist) %>%
  t() %>%
  as.data.frame() %>%
  dplyr::rename(functional = 1) %>%
  rownames_to_column(var = "sample") %>%
  mutate(functional = if_else(is.nan(functional), 1, functional))

# Merge all metrics
alpha_div <- richness %>%
  full_join(neutral, by = join_by(sample == sample)) %>%
  full_join(phylogenetic, by = join_by(sample == sample)) %>%
  full_join(functional, by = join_by(sample == sample))

6.1.1 Wild samples

alpha_div %>%
  pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
  left_join(., sample_metadata, by = join_by(sample == Tube_code)) %>%
  filter(time_point=="0_Wild") %>%
  mutate(metric=factor(metric,levels=c("richness","neutral","phylogenetic","functional"))) %>%
      ggplot(aes(y = value, x = Population, group=Population, color=Population, fill=Population)) +
      geom_boxplot(outlier.shape = NA) +
      geom_jitter(alpha=0.5) +
      scale_color_manual(name="Population",
          breaks=c("Cold_wet","Hot_dry"),
          labels=c("Cold","Hot"),
          values=c('#008080', "#d57d2c")) +
      scale_fill_manual(name="Population",
          breaks=c("Cold_wet","Hot_dry"),
          labels=c("Cold","Hot"),
          values=c('#00808050', "#d57d2c50")) +
      facet_wrap(. ~ metric,scales = "free") +
      coord_cartesian(xlim = c(1, NA)) +
      stat_compare_means(size=3, label.x=.58) +
      theme_classic() +
        theme(
    strip.background = element_blank(),
    panel.grid.minor.x = element_line(size = .1, color = "grey"),
    axis.title.x = element_blank(),
    axis.title.y = element_text(size=10),
    axis.text.x = element_text(angle = 45, hjust = 1),
    # Increase plot size
    plot.title = element_text(size = 10),
    axis.text = element_text(size = 8),
    axis.title = element_text(size = 8)
      ) +
  ylab("Alpha diversity")

6.1.2 Acclimation samples

alpha_div %>%
  pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
  left_join(., sample_metadata, by = join_by(sample == Tube_code)) %>%
  filter(time_point=="1_Acclimation") %>%
  mutate(metric=factor(metric,levels=c("richness","neutral","phylogenetic","functional"))) %>%
      ggplot(aes(y = value, x = Population, group=Population, color=Population, fill=Population)) +
      geom_boxplot(outlier.shape = NA) +
      geom_jitter(alpha=0.5) +
      scale_color_manual(name="Population",
          breaks=c("Cold_wet","Hot_dry"),
          labels=c("Cold","Hot"),
          values=c('#008080', "#d57d2c")) +
      scale_fill_manual(name="Population",
          breaks=c("Cold_wet","Hot_dry"),
          labels=c("Cold","Hot"),
          values=c('#00808050', "#d57d2c50")) +
      facet_wrap(. ~ metric,scales = "free") +
      coord_cartesian(xlim = c(1, NA)) +
      stat_compare_means(size=3, label.x=.58) +
      theme_classic() +
        theme(
    strip.background = element_blank(),
    panel.grid.minor.x = element_line(size = .1, color = "grey"),
    axis.title.x = element_blank(),
    axis.title.y = element_text(size=10),
    axis.text.x = element_text(angle = 45, hjust = 1),
    # Increase plot size
    plot.title = element_text(size = 10),
    axis.text = element_text(size = 8),
    axis.title = element_text(size = 8)
      ) +
  ylab("Alpha diversity")

6.1.3 Post-Transplant_1 samples

bxp1<-alpha_div %>%
  pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
  left_join(., sample_metadata, by = join_by(sample == Tube_code)) %>%
  filter(time_point=="5_Post-FMT1") %>%
  mutate(metric=factor(metric,levels=c("richness","neutral","phylogenetic","functional"))) %>%
      ggplot(aes(y = value, x = type, group=type, color=type, fill=type)) +
      geom_boxplot(outlier.shape = NA) +
      geom_jitter(alpha=0.5) +
      scale_color_manual(name="Type",
          breaks=c("Control","Hot_control", "Treatment"),
          labels=c("Cold-Cold","Hot-Hot", "Cold-Hot"),
          values=c("#4477AA","#d57d2c", "#76b183")) +
      scale_fill_manual(name="Type",
          breaks=c("Control","Hot_control", "Treatment"),
          labels=c("Cold-Cold","Hot-Hot", "Cold-Hot"),
          values=c("#4477AA50","#d57d2c50","#76b18350")) +
        scale_x_discrete(labels = c("Control" = "Cold-Cold", "Hot_control" = "Hot-Hot", "Treatment" = "Cold-Hot")) +
      facet_wrap(. ~ metric,scales = "free") +
      coord_cartesian(xlim = c(1, NA)) +
      stat_compare_means(size=3, label.x=.7) +
      theme_classic() +
        theme(
    strip.background = element_blank(),
    panel.grid.minor.x = element_line(size = .1, color = "grey"),
    axis.title.x = element_blank(),
    axis.title.y = element_text(size=10),
    axis.text.x = element_text(angle = 45, hjust = 1),
    # Increase plot size
    plot.title = element_text(size = 10),
    axis.text = element_text(size = 8),
    axis.title = element_text(size = 8)
      ) +
  ylab("Alpha diversity")


# Add significance bars with p-values
bxp1 + 
  geom_signif(comparisons = list(c("Control", "Hot_control"), c("Control", "Treatment"), c("Hot_control", "Treatment")),
              map_signif_level = TRUE)

6.1.3.1 Richness

alpha_div_meta <- alpha_div %>%
  left_join(sample_metadata, by = join_by(sample == Tube_code)) %>%
  filter(time_point =="5_Post-FMT1") 


Modelq0GLMMNB <- glmer.nb(richness ~ type+(1|individual), data = alpha_div_meta)
summary(Modelq0GLMMNB)
Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
 Family: Negative Binomial(15311.33)  ( log )
Formula: richness ~ type + (1 | individual)
   Data: alpha_div_meta

     AIC      BIC   logLik deviance df.resid 
   248.5    254.8   -119.3    238.5       21 

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-1.35323 -0.19624  0.00957  0.16293  0.31217 

Random effects:
 Groups     Name        Variance Std.Dev.
 individual (Intercept) 0.2757   0.5251  
Number of obs: 26, groups:  individual, 26

Fixed effects:
                Estimate Std. Error z value Pr(>|z|)    
(Intercept)       3.5033     0.1854  18.893   <2e-16 ***
typeHot_control   0.5046     0.2590   1.949   0.0513 .  
typeTreatment     0.2199     0.2692   0.817   0.4140    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
            (Intr) typHt_
typHt_cntrl -0.716       
typeTretmnt -0.686  0.491

6.1.3.2 Neutral

Model_neutral <- lme(fixed = neutral ~ type, data = alpha_div_meta,
               random = ~ 1 | individual)
summary(Model_neutral)
Linear mixed-effects model fit by REML
  Data: alpha_div_meta 
       AIC      BIC    logLik
  188.7488 194.4263 -89.37439

Random effects:
 Formula: ~1 | individual
        (Intercept) Residual
StdDev:    9.586725 3.595032

Fixed effects:  neutral ~ type 
                    Value Std.Error DF  t-value p-value
(Intercept)     17.359333  3.412877 23 5.086423  0.0000
typeHot_control  8.380488  4.826537 23 1.736336  0.0959
typeTreatment    7.504392  4.975080 23 1.508396  0.1451
 Correlation: 
                (Intr) typHt_
typeHot_control -0.707       
typeTreatment   -0.686  0.485

Standardized Within-Group Residuals:
        Min          Q1         Med          Q3         Max 
-0.70984520 -0.23549305 -0.03481478  0.24385684  0.54281137 

Number of Observations: 26
Number of Groups: 26 

6.1.3.3 Phylogenetic

Model_phylo <- lme(fixed = phylogenetic ~ Population, data = alpha_div_meta,
               random = ~ 1 | individual)
summary(Model_phylo)
Linear mixed-effects model fit by REML
  Data: alpha_div_meta 
       AIC      BIC    logLik
  72.19064 76.90286 -32.09532

Random effects:
 Formula: ~1 | individual
        (Intercept)  Residual
StdDev:   0.7770718 0.2914019

Fixed effects:  phylogenetic ~ Population 
                      Value Std.Error DF  t-value p-value
(Intercept)        4.307539 0.2012835 24 21.40036  0.0000
PopulationHot_dry -0.036754 0.3421161 24 -0.10743  0.9153
 Correlation: 
                  (Intr)
PopulationHot_dry -0.588

Standardized Within-Group Residuals:
        Min          Q1         Med          Q3         Max 
-0.71953431 -0.17146063  0.01286462  0.19655255  0.64008560 

Number of Observations: 26
Number of Groups: 26 

6.1.3.4 Functional

Model_func <- lme(fixed = functional ~ type, data = alpha_div_meta,
               random = ~ 1 | individual)
summary(Model_func)
Linear mixed-effects model fit by REML
  Data: alpha_div_meta 
        AIC       BIC   logLik
  -27.23436 -21.55689 18.61718

Random effects:
 Formula: ~1 | individual
        (Intercept)   Residual
StdDev:  0.08760598 0.03285224

Fixed effects:  functional ~ type 
                    Value  Std.Error DF  t-value p-value
(Intercept)     1.4854053 0.03118774 23 47.62786  0.0000
typeHot_control 0.0073281 0.04410613 23  0.16615  0.8695
typeTreatment   0.0032476 0.04546356 23  0.07143  0.9437
 Correlation: 
                (Intr) typHt_
typeHot_control -0.707       
typeTreatment   -0.686  0.485

Standardized Within-Group Residuals:
        Min          Q1         Med          Q3         Max 
-0.73243259 -0.23582828 -0.03001854  0.23261622  0.61672731 

Number of Observations: 26
Number of Groups: 26 

6.1.4 Post-Transplant_2 samples

bxp1<-alpha_div %>%
  pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
  left_join(., sample_metadata, by = join_by(sample == Tube_code)) %>%
  filter(time_point=="6_Post-FMT2") %>%
  mutate(metric=factor(metric,levels=c("richness","neutral","phylogenetic","functional"))) %>%
      ggplot(aes(y = value, x = type, group=type, color=type, fill=type)) +
      geom_boxplot(outlier.shape = NA) +
      geom_jitter(alpha=0.5) +
      scale_color_manual(name="Type",
          breaks=c("Control","Hot_control", "Treatment"),
          labels=c("Cold-control","Warm-control", "Cold-intervention"),
          values=c("#4477AA","#d57d2c", "#76b183")) +
      scale_fill_manual(name="Type",
          breaks=c("Control","Hot_control", "Treatment"),
          labels=c("Cold-control","Warm-control", "Cold-intervention"),
          values=c("#4477AA50","#d57d2c50","#76b18350")) +
          scale_x_discrete(labels = c("Control" = "Cold-Cold", "Hot_control" = "Hot-Hot", "Treatment" = "Cold-Hot")) +
      facet_wrap(. ~ metric,scales = "free") +
      coord_cartesian(xlim = c(1, NA)) +
      stat_compare_means(size=3, label.x=.7) +
      theme_classic() +
        theme(
    strip.background = element_blank(),
    panel.grid.minor.x = element_line(size = .1, color = "grey"),
    axis.title.x = element_blank(),
    axis.title.y = element_text(size=10),
    axis.text.x = element_text(angle = 45, hjust = 1),
    # Increase plot size
    plot.title = element_text(size = 10),
    axis.text = element_text(size = 8),
    axis.title = element_text(size = 8)
      ) +
  ylab("Alpha diversity")

# Add significance bars with p-values
bxp1 + 
  geom_signif(comparisons = list(c("Control", "Hot_control"), c("Control", "Treatment"), c("Hot_control", "Treatment")),
              map_signif_level = TRUE)

6.2 Beta diversity

beta_q0n <- genome_counts_filt %>%
  column_to_rownames(., "genome") %>%
  hillpair(., q = 0)

beta_q1n <- genome_counts_filt %>%
  column_to_rownames(., "genome") %>%
  hillpair(., q = 1)

beta_q1p <- genome_counts_filt %>%
  column_to_rownames(., "genome") %>%
  hillpair(., q = 1, tree = genome_tree)

beta_q1f <- genome_counts_filt %>%
  column_to_rownames(., "genome") %>%
  hillpair(., q = 1, dist = dist)

6.3 Permanovas

6.3.0.1 Load required data

6.3.1 1. Are the wild populations similar?

6.3.1.1 Wild: cold vs warm adapted lizards

6.3.1.1.1 Number of samples used
[1] 27

6.3.1.2 Richness

richness <- as.matrix(beta_q0n$S)
richness <- as.dist(richness[rownames(richness) %in% samples_to_keep,
               colnames(richness) %in% samples_to_keep])
betadisper(richness, subset_meta$Population) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df   Sum Sq  Mean Sq      F N.Perm Pr(>F)
Groups     1 0.000012 0.000012 0.0012    999  0.978
Residuals 25 0.257281 0.010291                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
         Cold_wet Hot_dry
Cold_wet             0.98
Hot_dry   0.97302        
adonis2(richness ~ Population,
        data = subset_meta %>% arrange(match(Tube_code,labels(richness))),
        permutations = 999) %>%
        broom::tidy() %>%
        tt()
term df SumOfSqs R2 statistic p.value
Model 1 1.542719 0.2095041 6.625717 0.001
Residual 25 5.820951 0.7904959 NA NA
Total 26 7.363669 1.0000000 NA NA

6.3.1.3 Neutral

neutral <- as.matrix(beta_q1n$S)
neutral <- as.dist(neutral[rownames(neutral) %in% samples_to_keep,
               colnames(neutral) %in% samples_to_keep])
betadisper(neutral, subset_meta$Population) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df   Sum Sq   Mean Sq      F N.Perm Pr(>F)
Groups     1 0.000048 0.0000476 0.0044    999  0.951
Residuals 25 0.270114 0.0108046                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
         Cold_wet Hot_dry
Cold_wet            0.954
Hot_dry   0.94763        
adonis2(neutral ~ Population,
        data = subset_meta %>% arrange(match(Tube_code,labels(neutral))),
        permutations = 999) %>%
        broom::tidy() %>%
        tt()
term df SumOfSqs R2 statistic p.value
Model 1 1.918266 0.2608511 8.822682 0.001
Residual 25 5.435610 0.7391489 NA NA
Total 26 7.353876 1.0000000 NA NA

6.3.1.4 Phylogenetic

phylo <- as.matrix(beta_q1p$S)
phylo <- as.dist(phylo[rownames(phylo) %in% samples_to_keep,
               colnames(phylo) %in% samples_to_keep])
betadisper(phylo, subset_meta$Population) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq  Mean Sq      F N.Perm Pr(>F)
Groups     1 0.03585 0.035847 2.4912    999  0.133
Residuals 25 0.35973 0.014389                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
         Cold_wet Hot_dry
Cold_wet            0.126
Hot_dry   0.12705        
adonis2(phylo ~ Population,
        data = subset_meta %>% arrange(match(Tube_code,labels(phylo))),
        permutations = 999) %>%
        broom::tidy() %>%
        tt()
term df SumOfSqs R2 statistic p.value
Model 1 0.3218613 0.2162815 6.899207 0.001
Residual 25 1.1662981 0.7837185 NA NA
Total 26 1.4881594 1.0000000 NA NA

6.3.1.5 Functional

func <- as.matrix(beta_q1f$S)
func <- as.dist(func[rownames(func) %in% samples_to_keep,
               colnames(func) %in% samples_to_keep])
betadisper(func, subset_meta$Population) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df   Sum Sq  Mean Sq      F N.Perm Pr(>F)
Groups     1 0.018367 0.018367 1.5597    999  0.209
Residuals 25 0.294402 0.011776                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
         Cold_wet Hot_dry
Cold_wet            0.224
Hot_dry   0.22328        
adonis2(func ~ Population,
        data = subset_meta %>% arrange(match(Tube_code,labels(func))),
        permutations = 999) %>%
        broom::tidy() %>%
        tt()
term df SumOfSqs R2 statistic p.value
Model 1 0.08585779 0.172879 5.225323 0.063
Residual 25 0.41077745 0.827121 NA NA
Total 26 0.49663525 1.000000 NA NA
beta_q0n_nmds_wild <- richness %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE, trace=FALSE) %>%
                vegan::scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(subset_meta, by = join_by(sample == Tube_code))

beta_q1n_nmds_wild <- neutral %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE, trace=FALSE) %>%
                vegan::scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(subset_meta, by = join_by(sample == Tube_code))

beta_q1p_nmds_wild <- phylo %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                vegan::scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(subset_meta, by = join_by(sample == Tube_code))

beta_q1f_nmds_wild <- func %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                vegan::scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(subset_meta, by = join_by(sample == Tube_code))

6.3.2 2. Effect of acclimation

6.3.2.0.1 Number of samples used
[1] 26

6.3.2.1 Richness

richness_accli <- as.matrix(beta_q0n$S)
richness_accli <- as.dist(richness_accli[rownames(richness_accli) %in% samples_to_keep_accli,
               colnames(richness_accli) %in% samples_to_keep_accli])
betadisper(richness_accli, subset_meta_accli$Population) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df   Sum Sq  Mean Sq      F N.Perm Pr(>F)   
Groups     1 0.093187 0.093187 11.812    999  0.004 **
Residuals 24 0.189340 0.007889                        
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
          Cold_wet Hot_dry
Cold_wet             0.002
Hot_dry  0.0021532        
adonis2(richness_accli ~ Population,
        data = subset_meta_accli %>% arrange(match(Tube_code,labels(richness_accli))),
        permutations = 999) %>%
        broom::tidy() %>%
        tt()
term df SumOfSqs R2 statistic p.value
Model 1 1.639630 0.1929088 5.736415 0.001
Residual 24 6.859879 0.8070912 NA NA
Total 25 8.499509 1.0000000 NA NA

6.3.2.2 Neutral

neutral_accli <- as.matrix(beta_q1n$S)
neutral_accli <- as.dist(neutral_accli[rownames(neutral_accli) %in% samples_to_keep_accli,
               colnames(neutral_accli) %in% samples_to_keep_accli])
betadisper(neutral_accli, subset_meta_accli$Population) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq  Mean Sq      F N.Perm Pr(>F)  
Groups     1 0.05603 0.056026 4.1918    999  0.041 *
Residuals 24 0.32077 0.013365                       
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
         Cold_wet Hot_dry
Cold_wet            0.043
Hot_dry  0.051717        
adonis2(neutral_accli ~ Population,
        data = subset_meta_accli %>% arrange(match(Tube_code,labels(neutral_accli))),
        permutations = 999) %>%
        broom::tidy() %>%
        tt()
term df SumOfSqs R2 statistic p.value
Model 1 1.961192 0.2493171 7.970889 0.001
Residual 24 5.905063 0.7506829 NA NA
Total 25 7.866255 1.0000000 NA NA

6.3.2.3 Phylogenetic

phylo_accli <- as.matrix(beta_q1p$S)
phylo_accli <- as.dist(phylo_accli[rownames(phylo_accli) %in% samples_to_keep_accli,
               colnames(phylo_accli) %in% samples_to_keep_accli])
betadisper(phylo_accli, subset_meta_accli$Population) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq  Mean Sq      F N.Perm Pr(>F)
Groups     1 0.03637 0.036365 2.3087    999  0.142
Residuals 24 0.37804 0.015752                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
         Cold_wet Hot_dry
Cold_wet             0.15
Hot_dry   0.14172        
adonis2(phylo_accli ~ Population,
        data = subset_meta_accli %>% arrange(match(Tube_code,labels(phylo_accli))),
        permutations = 999) %>%
        broom::tidy() %>%
        tt()
term df SumOfSqs R2 statistic p.value
Model 1 0.2113395 0.1379515 3.84066 0.009
Residual 24 1.3206449 0.8620485 NA NA
Total 25 1.5319844 1.0000000 NA NA

6.3.2.4 Functional

func_accli <- as.matrix(beta_q1f$S)
func_accli <- as.dist(func_accli[rownames(func_accli) %in% samples_to_keep_accli,
               colnames(func_accli) %in% samples_to_keep_accli])
betadisper(func_accli, subset_meta_accli$Population) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq   Mean Sq      F N.Perm Pr(>F)
Groups     1 0.00471 0.0047108 0.2042    999  0.662
Residuals 24 0.55370 0.0230708                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
         Cold_wet Hot_dry
Cold_wet            0.664
Hot_dry   0.65542        
adonis2(func_accli ~ Population,
        data = subset_meta_accli %>% arrange(match(Tube_code,labels(func_accli))),
        permutations = 999) %>%
        broom::tidy() %>%
        tt()
term df SumOfSqs R2 statistic p.value
Model 1 0.00470883 0.006952167 0.1680201 0.614
Residual 24 0.67260942 0.993047833 NA NA
Total 25 0.67731825 1.000000000 NA NA
beta_q0n_nmds_accli <- richness_accli %>%
  metaMDS(.,trymax = 500, k=2, verbosity=FALSE, trace=FALSE) %>%
  vegan::scores() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(subset_meta_accli, by = join_by(sample == Tube_code))

beta_q1n_nmds_accli <- neutral_accli %>%
  metaMDS(.,trymax = 500, k=2, verbosity=FALSE, trace=FALSE) %>%
  vegan::scores() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(subset_meta_accli, by = join_by(sample == Tube_code))

beta_q1p_nmds_accli <- phylo_accli %>%
  metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
  vegan::scores() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(subset_meta_accli, by = join_by(sample == Tube_code))

beta_q1f_nmds_accli <- func_accli %>%
  metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
  vegan::scores() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(subset_meta_accli, by = join_by(sample == Tube_code))

6.3.3 3. Comparison between Wild and Acclimation

6.3.3.0.1 Number of samples used
[1] 53
6.3.3.0.2 Richness
richness_accli1 <- as.matrix(beta_q0n$S)
richness_accli1 <- as.dist(richness_accli1[rownames(richness_accli1) %in% samples_to_keep_accli1,
               colnames(richness_accli1) %in% samples_to_keep_accli1])
betadisper(richness_accli1, subset_meta_accli1$Population) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq  Mean Sq      F N.Perm Pr(>F)  
Groups     1 0.03286 0.032865 2.9698    999  0.094 .
Residuals 51 0.56438 0.011066                       
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
         Cold_wet Hot_dry
Cold_wet            0.083
Hot_dry   0.09089        
adonis2(richness_accli1 ~ Population*time_point,
        data = subset_meta_accli1 %>% arrange(match(Tube_code,labels(richness_accli1))),
        permutations = 999,
        strata = subset_meta_accli1 %>% arrange(match(Tube_code,labels(richness_accli1))) %>% pull(individual)) %>%
        broom::tidy() %>%
        tt()
term df SumOfSqs R2 statistic p.value
Model 3 3.762933 0.2288365 4.846784 0.001
Residual 49 12.680830 0.7711635 NA NA
Total 52 16.443763 1.0000000 NA NA
#Arrange of metadata dataframe
subset_meta_accli1_arrange <- column_to_rownames(subset_meta_accli1, "Tube_code")
subset_meta_accli1_arrange<-subset_meta_accli1_arrange[labels(richness_accli1),]

pairwise<-pairwise.adonis(richness_accli1,subset_meta_accli1_arrange$Population_time, perm=999)
pairwise%>%
  tt()
pairs Df SumsOfSqs F.Model R2 p.value p.adjusted sig
Cold_wet.1_Acclimation vs Hot_dry.1_Acclimation 1 1.6396299 5.736415 0.19290877 0.001 0.006 *
Cold_wet.1_Acclimation vs Cold_wet.0_Wild 1 0.5273732 1.906862 0.05462715 0.010 0.060
Cold_wet.1_Acclimation vs Hot_dry.0_Wild 1 1.5558412 5.190707 0.17782052 0.001 0.006 *
Hot_dry.1_Acclimation vs Cold_wet.0_Wild 1 1.8259388 8.319131 0.24968031 0.001 0.006 *
Hot_dry.1_Acclimation vs Hot_dry.0_Wild 1 0.3856333 1.736034 0.09788177 0.009 0.054
Cold_wet.0_Wild vs Hot_dry.0_Wild 1 1.5427188 6.625717 0.20950408 0.001 0.006 *
6.3.3.0.3 Neutral
neutral_accli1 <- as.matrix(beta_q1n$S)
neutral_accli1 <- as.dist(neutral_accli1[rownames(neutral_accli1) %in% samples_to_keep_accli1,
               colnames(neutral_accli1) %in% samples_to_keep_accli1])
betadisper(neutral_accli1, subset_meta_accli1$Population) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq  Mean Sq      F N.Perm Pr(>F)
Groups     1 0.01431 0.014310 1.1177    999  0.313
Residuals 51 0.65296 0.012803                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
         Cold_wet Hot_dry
Cold_wet            0.297
Hot_dry   0.29539        
adonis2(neutral_accli1 ~ Population*time_point,
        data = subset_meta_accli1 %>% arrange(match(Tube_code,labels(neutral_accli1))),
        permutations = 999,
        strata = subset_meta_accli1 %>% arrange(match(Tube_code,labels(neutral_accli1))) %>% pull(individual)) %>%
        broom::tidy() %>%
        tt()
term df SumOfSqs R2 statistic p.value
Model 3 4.735485 0.2945657 6.820252 0.001
Residual 49 11.340673 0.7054343 NA NA
Total 52 16.076158 1.0000000 NA NA
pairwise<-pairwise.adonis(neutral_accli1,subset_meta_accli1_arrange$Population_time, perm=999)
pairwise%>%
  tt()
pairs Df SumsOfSqs F.Model R2 p.value p.adjusted sig
Cold_wet.1_Acclimation vs Hot_dry.1_Acclimation 1 1.9611917 7.970889 0.24931708 0.001 0.006 *
Cold_wet.1_Acclimation vs Cold_wet.0_Wild 1 0.7466371 3.064914 0.08498327 0.001 0.006 *
Cold_wet.1_Acclimation vs Hot_dry.0_Wild 1 2.1421077 8.176866 0.25412251 0.001 0.006 *
Hot_dry.1_Acclimation vs Cold_wet.0_Wild 1 2.0371666 10.078295 0.28730857 0.001 0.006 *
Hot_dry.1_Acclimation vs Hot_dry.0_Wild 1 0.6314393 3.060027 0.16054681 0.001 0.006 *
Cold_wet.0_Wild vs Hot_dry.0_Wild 1 1.9182663 8.822682 0.26085105 0.001 0.006 *
6.3.3.0.4 Phylogenetic
phylo_accli1 <- as.matrix(beta_q1p$S)
phylo_accli1 <- as.dist(phylo_accli1[rownames(phylo_accli1) %in% samples_to_keep_accli1,
               colnames(phylo_accli1) %in% samples_to_keep_accli1])
betadisper(phylo_accli1, subset_meta_accli1$Population) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq   Mean Sq     F N.Perm Pr(>F)
Groups     1 0.00001 0.0000111 6e-04    999  0.984
Residuals 51 0.89017 0.0174543                    

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
         Cold_wet Hot_dry
Cold_wet            0.985
Hot_dry   0.97998        
adonis2(phylo_accli1 ~ Population*time_point,
        data = subset_meta_accli1 %>% arrange(match(Tube_code,labels(phylo_accli1))),
        permutations = 999,
        strata = subset_meta_accli1 %>% arrange(match(Tube_code,labels(phylo_accli1))) %>% pull(individual)) %>%
        broom::tidy() %>%
        tt()
term df SumOfSqs R2 statistic p.value
Model 3 0.7622566 0.2345983 5.006223 0.001
Residual 49 2.4869430 0.7654017 NA NA
Total 52 3.2491997 1.0000000 NA NA
pairwise<-pairwise.adonis(phylo_accli1,subset_meta_accli1_arrange$Population_time, perm=999)
pairwise%>%
  tt()
pairs Df SumsOfSqs F.Model R2 p.value p.adjusted sig
Cold_wet.1_Acclimation vs Hot_dry.1_Acclimation 1 0.2113395 3.840660 0.1379515 0.010 0.060
Cold_wet.1_Acclimation vs Cold_wet.0_Wild 1 0.2076830 4.216762 0.1133028 0.003 0.018 .
Cold_wet.1_Acclimation vs Hot_dry.0_Wild 1 0.3635093 4.799786 0.1666605 0.001 0.006 *
Hot_dry.1_Acclimation vs Cold_wet.0_Wild 1 0.2121479 7.924059 0.2406769 0.002 0.012 .
Hot_dry.1_Acclimation vs Hot_dry.0_Wild 1 0.2092433 3.885515 0.1953942 0.001 0.006 *
Cold_wet.0_Wild vs Hot_dry.0_Wild 1 0.3218613 6.899207 0.2162815 0.001 0.006 *
6.3.3.0.5 Functional
func_accli1 <- as.matrix(beta_q1f$S)
func_accli1 <- as.dist(func_accli1[rownames(func_accli1) %in% samples_to_keep_accli1,
               colnames(func_accli1) %in% samples_to_keep_accli1])
betadisper(func_accli1, subset_meta_accli1$Population) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq  Mean Sq      F N.Perm Pr(>F)
Groups     1 0.01354 0.013538 0.7762    999   0.41
Residuals 51 0.88953 0.017442                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
         Cold_wet Hot_dry
Cold_wet            0.408
Hot_dry   0.38245        
adonis2(func_accli1 ~ Population*time_point,
        data = subset_meta_accli1 %>% arrange(match(Tube_code,labels(func_accli1))),
        permutations = 999,
        strata = subset_meta_accli1 %>% arrange(match(Tube_code,labels(func_accli1))) %>% pull(individual)) %>%
        broom::tidy() %>%
        tt()
term df SumOfSqs R2 statistic p.value
Model 3 0.1059556 0.08908753 1.597405 0.206
Residual 49 1.0833869 0.91091247 NA NA
Total 52 1.1893425 1.00000000 NA NA
pairwise<-pairwise.adonis(func_accli1,subset_meta_accli1_arrange$Population_time, perm=999)
pairwise%>%
  tt()
pairs Df SumsOfSqs F.Model R2 p.value p.adjusted sig
Cold_wet.1_Acclimation vs Hot_dry.1_Acclimation 1 0.004708830 0.16802012 0.006952167 0.603 1.000
Cold_wet.1_Acclimation vs Cold_wet.0_Wild 1 0.001001660 0.03572249 0.001081329 0.691 1.000
Cold_wet.1_Acclimation vs Hot_dry.0_Wild 1 0.088628476 3.78330596 0.136171914 0.092 0.552
Hot_dry.1_Acclimation vs Cold_wet.0_Wild 1 0.002762199 0.13250295 0.005272175 0.605 1.000
Hot_dry.1_Acclimation vs Hot_dry.0_Wild 1 0.042282101 4.27996390 0.211043961 0.096 0.576
Cold_wet.0_Wild vs Hot_dry.0_Wild 1 0.085857794 5.22532292 0.172878978 0.055 0.330
beta_richness_nmds_accli1 <- richness_accli1 %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(subset_meta_accli1, by = c("sample" = "Tube_code"))

beta_neutral_nmds_accli1 <- neutral_accli1 %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(subset_meta_accli1, by = c("sample" = "Tube_code"))

beta_phylo_nmds_accli1 <- phylo_accli1 %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(subset_meta_accli1, by = join_by(sample == Tube_code))

beta_func_nmds_accli1 <- func_accli1 %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(subset_meta_accli1, by = join_by(sample == Tube_code))

6.3.4 4. Effect of FMT on microbiota community

6.3.4.1 Comparison between Acclimation vs Post-FMT1

6.3.4.1.1 Number of samples used
[1] 52
6.3.4.1.2 Richness
richness_post3 <- as.matrix(beta_q0n$S)
richness_post3 <- as.dist(richness_post3[rownames(richness_post3) %in% samples_to_keep_post3,
               colnames(richness_post3) %in% samples_to_keep_post3])
betadisper(richness_post3, subset_meta_post3$Population) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq  Mean Sq      F N.Perm Pr(>F)    
Groups     1 0.10843 0.108427 25.578    999  0.001 ***
Residuals 50 0.21195 0.004239                         
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
           Cold_wet Hot_dry
Cold_wet              0.001
Hot_dry  6.0928e-06        
adonis2(richness_post3 ~ Population*time_point,
        data = subset_meta_post3 %>% arrange(match(Tube_code,labels(richness_post3))),
        permutations = 999,
        strata = subset_meta_post3 %>% arrange(match(Tube_code,labels(richness_post3))) %>% pull(individual)) %>%
        broom::tidy() %>%
        tt()
term df SumOfSqs R2 statistic p.value
Model 3 3.472191 0.193401 3.836375 0.001
Residual 48 14.481131 0.806599 NA NA
Total 51 17.953321 1.000000 NA NA
#Arrange of metadata dataframe
subset_meta_post3_arrange <- column_to_rownames(subset_meta_post3, "Tube_code")
subset_meta_post3_arrange<-subset_meta_post3_arrange[labels(richness_post3),]

pairwise <- pairwise.adonis(richness_post3, subset_meta_post3_arrange$type_time, perm=999)
pairwise%>%
  tt()
pairs Df SumsOfSqs F.Model R2 p.value p.adjusted sig
Control.1_Acclimation vs Treatment.1_Acclimation 1 0.3657243 1.123239 0.06966584 0.258 1.000
Control.1_Acclimation vs Hot_control.1_Acclimation 1 1.2605860 4.943717 0.24788342 0.001 0.015 .
Control.1_Acclimation vs Control.5_Post-FMT1 1 0.6804988 2.100611 0.12283837 0.001 0.015 .
Control.1_Acclimation vs Treatment.5_Post-FMT1 1 0.8481258 2.673290 0.16033371 0.001 0.015 .
Control.1_Acclimation vs Hot_control.5_Post-FMT1 1 1.1180661 3.809275 0.20252111 0.001 0.015 .
Treatment.1_Acclimation vs Hot_control.1_Acclimation 1 1.3606630 5.087152 0.24124415 0.001 0.015 .
Treatment.1_Acclimation vs Control.5_Post-FMT1 1 0.7216200 2.172734 0.11956009 0.001 0.015 .
Treatment.1_Acclimation vs Treatment.5_Post-FMT1 1 0.9551308 2.926054 0.16322910 0.001 0.015 .
Treatment.1_Acclimation vs Hot_control.5_Post-FMT1 1 1.2263345 4.039487 0.20157637 0.001 0.015 .
Hot_control.1_Acclimation vs Control.5_Post-FMT1 1 1.4319792 5.384836 0.25180628 0.001 0.015 .
Hot_control.1_Acclimation vs Treatment.5_Post-FMT1 1 0.8172413 3.194690 0.17558364 0.001 0.015 .
Hot_control.1_Acclimation vs Hot_control.5_Post-FMT1 1 0.5796135 2.441615 0.13239702 0.001 0.015 .
Control.5_Post-FMT1 vs Treatment.5_Post-FMT1 1 0.5615418 1.729004 0.10335366 0.013 0.195
Control.5_Post-FMT1 vs Hot_control.5_Post-FMT1 1 0.8438429 2.793772 0.14865413 0.001 0.015 .
Treatment.5_Post-FMT1 vs Hot_control.5_Post-FMT1 1 0.3734921 1.268929 0.07799710 0.108 1.000
6.3.4.1.3 Neutral
neutral_post3 <- as.matrix(beta_q1n$S)
neutral_post3 <- as.dist(neutral_post3[rownames(neutral_post3) %in% samples_to_keep_post3,
               colnames(neutral_post3) %in% samples_to_keep_post3])
betadisper(neutral_post3, subset_meta_post3$Population) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq  Mean Sq      F N.Perm Pr(>F)   
Groups     1 0.09089 0.090889 12.898    999  0.002 **
Residuals 50 0.35233 0.007047                        
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
           Cold_wet Hot_dry
Cold_wet              0.001
Hot_dry  0.00074926        
adonis2(neutral_post3 ~ Population*time_point,
        data = subset_meta_post3 %>% arrange(match(Tube_code,labels(neutral_post3))),
        permutations = 999,
        strata = subset_meta_post3 %>% arrange(match(Tube_code,labels(neutral_post3))) %>% pull(individual)) %>%
        broom::tidy() %>%
        tt()
term df SumOfSqs R2 statistic p.value
Model 3 4.443476 0.2628643 5.705637 0.001
Residual 48 12.460591 0.7371357 NA NA
Total 51 16.904067 1.0000000 NA NA
pairwise <- pairwise.adonis(neutral_post3, subset_meta_post3_arrange$type_time, perm=999)
pairwise%>%
  tt()
pairs Df SumsOfSqs F.Model R2 p.value p.adjusted sig
Control.1_Acclimation vs Treatment.1_Acclimation 1 0.2759858 0.9928976 0.06208366 0.457 1.000
Control.1_Acclimation vs Hot_control.1_Acclimation 1 1.4146464 6.5078325 0.30257965 0.001 0.015 .
Control.1_Acclimation vs Control.5_Post-FMT1 1 0.8153894 3.0970603 0.17113610 0.001 0.015 .
Control.1_Acclimation vs Treatment.5_Post-FMT1 1 1.1809241 4.4856470 0.24265567 0.001 0.015 .
Control.1_Acclimation vs Hot_control.5_Post-FMT1 1 1.4321524 5.7774260 0.27806264 0.001 0.015 .
Treatment.1_Acclimation vs Hot_control.1_Acclimation 1 1.6347704 6.8326887 0.29925029 0.001 0.015 .
Treatment.1_Acclimation vs Control.5_Post-FMT1 1 0.9517634 3.3715700 0.17404733 0.001 0.015 .
Treatment.1_Acclimation vs Treatment.5_Post-FMT1 1 1.3127773 4.6298256 0.23585668 0.002 0.030 .
Treatment.1_Acclimation vs Hot_control.5_Post-FMT1 1 1.6713369 6.2395460 0.28056085 0.001 0.015 .
Hot_control.1_Acclimation vs Control.5_Post-FMT1 1 1.5409781 6.8338056 0.29928456 0.001 0.015 .
Hot_control.1_Acclimation vs Treatment.5_Post-FMT1 1 0.9133614 4.0964534 0.21451383 0.001 0.015 .
Hot_control.1_Acclimation vs Hot_control.5_Post-FMT1 1 0.6954835 3.2951234 0.17077493 0.001 0.015 .
Control.5_Post-FMT1 vs Treatment.5_Post-FMT1 1 0.6051778 2.2508491 0.13047758 0.011 0.165
Control.5_Post-FMT1 vs Hot_control.5_Post-FMT1 1 1.0528902 4.1436369 0.20570451 0.002 0.030 .
Treatment.5_Post-FMT1 vs Hot_control.5_Post-FMT1 1 0.4150076 1.6372683 0.09840968 0.049 0.735
6.3.4.1.4 Phylogenetic
phylo_post3 <- as.matrix(beta_q1p$S)
phylo_post3 <- as.dist(phylo_post3[rownames(phylo_post3) %in% samples_to_keep_post3,
               colnames(phylo_post3) %in% samples_to_keep_post3])
betadisper(phylo_post3, subset_meta_post3$Population) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq  Mean Sq      F N.Perm Pr(>F)  
Groups     1 0.04088 0.040882 2.9254    999  0.097 .
Residuals 50 0.69874 0.013975                       
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
         Cold_wet Hot_dry
Cold_wet            0.099
Hot_dry  0.093394        
adonis2(phylo_post3 ~ Population*time_point,
        data = subset_meta_post3 %>% arrange(match(Tube_code,labels(phylo_post3))),
        permutations = 999,
        strata = subset_meta_post3 %>% arrange(match(Tube_code,labels(phylo_post3))) %>% pull(individual)) %>%
        broom::tidy() %>%
        tt()
term df SumOfSqs R2 statistic p.value
Model 3 0.6468577 0.2180259 4.461035 0.001
Residual 48 2.3200272 0.7819741 NA NA
Total 51 2.9668849 1.0000000 NA NA
pairwise <- pairwise.adonis(phylo_post3, subset_meta_post3_arrange$type_time, perm=999)
pairwise%>%
  tt()
pairs Df SumsOfSqs F.Model R2 p.value p.adjusted sig
Control.1_Acclimation vs Treatment.1_Acclimation 1 0.05690611 0.7893300 0.04999136 0.536 1.000
Control.1_Acclimation vs Hot_control.1_Acclimation 1 0.10866997 3.0547314 0.16919285 0.012 0.180
Control.1_Acclimation vs Control.5_Post-FMT1 1 0.11760287 2.9988032 0.16661126 0.015 0.225
Control.1_Acclimation vs Treatment.5_Post-FMT1 1 0.09624422 1.7247115 0.10968160 0.175 1.000
Control.1_Acclimation vs Hot_control.5_Post-FMT1 1 0.18586788 4.1904227 0.21836011 0.005 0.075
Treatment.1_Acclimation vs Hot_control.1_Acclimation 1 0.23108846 4.0521838 0.20208192 0.004 0.060
Treatment.1_Acclimation vs Control.5_Post-FMT1 1 0.26358465 4.3608960 0.21417997 0.004 0.060
Treatment.1_Acclimation vs Treatment.5_Post-FMT1 1 0.25319427 3.2738422 0.17915456 0.038 0.570
Treatment.1_Acclimation vs Hot_control.5_Post-FMT1 1 0.39050120 5.9837393 0.27218933 0.001 0.015 .
Hot_control.1_Acclimation vs Control.5_Post-FMT1 1 0.14203376 5.4200212 0.25303529 0.001 0.015 .
Hot_control.1_Acclimation vs Treatment.5_Post-FMT1 1 0.09666753 2.3682173 0.13635351 0.023 0.345
Hot_control.1_Acclimation vs Hot_control.5_Post-FMT1 1 0.09252600 2.9824958 0.15711821 0.011 0.165
Control.5_Post-FMT1 vs Treatment.5_Post-FMT1 1 0.01842535 0.4144162 0.02688498 0.737 1.000
Control.5_Post-FMT1 vs Hot_control.5_Post-FMT1 1 0.05987967 1.7387847 0.09802164 0.115 1.000
Treatment.5_Post-FMT1 vs Hot_control.5_Post-FMT1 1 0.03212966 0.6477782 0.04139746 0.694 1.000
6.3.4.1.5 Functional
func_post3 <- as.matrix(beta_q1f$S)
func_post3 <- as.dist(func_post3[rownames(func_post3) %in% samples_to_keep_post3,
               colnames(func_post3) %in% samples_to_keep_post3])
betadisper(func_post3, subset_meta_post3$Population) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq   Mean Sq      F N.Perm Pr(>F)
Groups     1 0.00125 0.0012477 0.0665    999  0.798
Residuals 50 0.93781 0.0187562                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
         Cold_wet Hot_dry
Cold_wet            0.799
Hot_dry   0.79753        
adonis2(func_post3 ~ Population*time_point,
        data = subset_meta_post3 %>% arrange(match(Tube_code,labels(func_post3))),
        permutations = 999,
        strata = subset_meta_post3 %>% arrange(match(Tube_code,labels(func_post3))) %>% pull(individual)) %>%
        broom::tidy() %>%
        tt()
term df SumOfSqs R2 statistic p.value
Model 3 0.03282552 0.02249363 0.3681798 0.363
Residual 48 1.42649950 0.97750637 NA NA
Total 51 1.45932502 1.00000000 NA NA
pairwise <- pairwise.adonis(func_post3, subset_meta_post3_arrange$type_time, perm=999)
pairwise%>%
  tt()
pairs Df SumsOfSqs F.Model R2 p.value p.adjusted sig
Control.1_Acclimation vs Treatment.1_Acclimation 1 0.022138495 0.64325180 0.041120082 0.406 1.00
Control.1_Acclimation vs Hot_control.1_Acclimation 1 0.017467416 0.54472723 0.035042572 0.440 1.00
Control.1_Acclimation vs Control.5_Post-FMT1 1 -0.007536319 -0.20662255 -0.013967233 0.906 1.00
Control.1_Acclimation vs Treatment.5_Post-FMT1 1 0.128319067 3.48711706 0.199410632 0.068 1.00
Control.1_Acclimation vs Hot_control.5_Post-FMT1 1 0.051950542 1.26548063 0.077801613 0.293 1.00
Treatment.1_Acclimation vs Hot_control.1_Acclimation 1 0.001373886 0.07238159 0.004503476 0.640 1.00
Treatment.1_Acclimation vs Control.5_Post-FMT1 1 0.002173211 0.09402475 0.005842215 0.632 1.00
Treatment.1_Acclimation vs Treatment.5_Post-FMT1 1 0.059119657 2.62461793 0.148917721 0.145 1.00
Treatment.1_Acclimation vs Hot_control.5_Post-FMT1 1 0.010911935 0.39816986 0.024281360 0.484 1.00
Hot_control.1_Acclimation vs Control.5_Post-FMT1 1 -0.002303582 -0.11016709 -0.006933181 0.748 1.00
Hot_control.1_Acclimation vs Treatment.5_Post-FMT1 1 0.058908140 2.91987617 0.162940644 0.165 1.00
Hot_control.1_Acclimation vs Hot_control.5_Post-FMT1 1 0.005904278 0.23427876 0.014431116 0.553 1.00
Control.5_Post-FMT1 vs Treatment.5_Post-FMT1 1 0.116146557 4.72479132 0.239535681 0.066 0.99
Control.5_Post-FMT1 vs Hot_control.5_Post-FMT1 1 0.050009298 1.70482607 0.096291602 0.237 1.00
Treatment.5_Post-FMT1 vs Hot_control.5_Post-FMT1 1 0.012358590 0.42381202 0.027477774 0.482 1.00
beta_richness_nmds_post3 <- richness_post3 %>%
  metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
  scores() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(subset_meta_post3, by = c("sample" = "Tube_code"))

beta_neutral_nmds_post3 <- neutral_post3 %>%
  metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
  scores() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(subset_meta_post3, by = c("sample" = "Tube_code"))

beta_phylo_nmds_post3 <- phylo_post3 %>%
  metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
  scores() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(subset_meta_post3, by = join_by(sample == Tube_code))

beta_func_nmds_post3 <- func_post3 %>%
  metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
  scores() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(subset_meta_post3, by = join_by(sample == Tube_code))

6.3.4.2 Comparison between Acclimation vs Post-FMT2

6.3.4.2.1 Number of samples used
[1] 53
6.3.4.2.2 Richness
richness_post4 <- as.matrix(beta_q0n$S)
richness_post4 <- as.dist(richness_post4[rownames(richness_post4) %in% samples_to_keep_post4,
               colnames(richness_post4) %in% samples_to_keep_post4])
betadisper(richness_post4, subset_meta_post4$Population) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq  Mean Sq      F N.Perm Pr(>F)    
Groups     1 0.07832 0.078322 11.371    999  0.001 ***
Residuals 51 0.35129 0.006888                         
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
          Cold_wet Hot_dry
Cold_wet             0.001
Hot_dry  0.0014303        
adonis2(richness_post4 ~ Population*time_point,
        data = subset_meta_post4 %>% arrange(match(Tube_code,labels(richness_post4))),
        permutations = 999,
        strata = subset_meta_post4 %>% arrange(match(Tube_code,labels(richness_post4))) %>% pull(individual)) %>%
        broom::tidy() %>%
        tt()
term df SumOfSqs R2 statistic p.value
Model 3 3.279021 0.1937289 3.924535 0.001
Residual 49 13.646799 0.8062711 NA NA
Total 52 16.925820 1.0000000 NA NA
#Arrange of metadata dataframe
subset_meta_post4_arrange <- column_to_rownames(subset_meta_post4, "Tube_code")
subset_meta_post4_arrange<-subset_meta_post4_arrange[labels(richness_post4),]

pairwise <- pairwise.adonis(richness_post4, subset_meta_post4_arrange$type_time, perm=999)
pairwise%>%
  tt()
pairs Df SumsOfSqs F.Model R2 p.value p.adjusted sig
Control.1_Acclimation vs Treatment.1_Acclimation 1 0.3657243 1.123239 0.06966584 0.236 1.000
Control.1_Acclimation vs Hot_control.1_Acclimation 1 1.2605860 4.943717 0.24788342 0.001 0.015 .
Control.1_Acclimation vs Treatment.6_Post-FMT2 1 0.8072604 2.940901 0.16392161 0.001 0.015 .
Control.1_Acclimation vs Control.6_Post-FMT2 1 0.4817387 1.660775 0.09968176 0.023 0.345
Control.1_Acclimation vs Hot_control.6_Post-FMT2 1 1.1179704 3.885459 0.20573812 0.001 0.015 .
Treatment.1_Acclimation vs Hot_control.1_Acclimation 1 1.3606630 5.087152 0.24124415 0.001 0.015 .
Treatment.1_Acclimation vs Treatment.6_Post-FMT2 1 0.9130048 3.195028 0.16645080 0.001 0.015 .
Treatment.1_Acclimation vs Control.6_Post-FMT2 1 0.5959230 1.984036 0.11032208 0.002 0.030 .
Treatment.1_Acclimation vs Hot_control.6_Post-FMT2 1 1.2747787 4.275366 0.21086503 0.001 0.015 .
Hot_control.1_Acclimation vs Treatment.6_Post-FMT2 1 0.6397330 2.913695 0.15405213 0.001 0.015 .
Hot_control.1_Acclimation vs Control.6_Post-FMT2 1 1.4575447 6.224524 0.28007456 0.001 0.015 .
Hot_control.1_Acclimation vs Hot_control.6_Post-FMT2 1 0.3276169 1.412318 0.08111028 0.037 0.555
Treatment.6_Post-FMT2 vs Control.6_Post-FMT2 1 0.6463814 2.560441 0.13795154 0.001 0.015 .
Treatment.6_Post-FMT2 vs Hot_control.6_Post-FMT2 1 0.4796256 1.916520 0.10696943 0.003 0.045 .
Control.6_Post-FMT2 vs Hot_control.6_Post-FMT2 1 1.1305044 4.268317 0.21059061 0.001 0.015 .
6.3.4.2.3 Neutral
neutral_post4 <- as.matrix(beta_q1n$S)
neutral_post4 <- as.dist(neutral_post4[rownames(neutral_post4) %in% samples_to_keep_post4,
               colnames(neutral_post4) %in% samples_to_keep_post4])
betadisper(neutral_post4, subset_meta_post4$Population) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq  Mean Sq      F N.Perm Pr(>F)   
Groups     1 0.07811 0.078108 9.6342    999  0.003 **
Residuals 51 0.41347 0.008107                        
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
          Cold_wet Hot_dry
Cold_wet             0.004
Hot_dry  0.0031133        
adonis2(neutral_post4 ~ Population*time_point,
        data = subset_meta_post4 %>% arrange(match(Tube_code,labels(neutral_post4))),
        permutations = 999,
        strata = subset_meta_post4 %>% arrange(match(Tube_code,labels(neutral_post4))) %>% pull(individual)) %>%
        broom::tidy() %>%
        tt()
term df SumOfSqs R2 statistic p.value
Model 3 3.823089 0.2357666 5.038847 0.001
Residual 49 12.392477 0.7642334 NA NA
Total 52 16.215567 1.0000000 NA NA
pairwise <- pairwise.adonis(neutral_post4, subset_meta_post4_arrange$type_time, perm=999)
pairwise%>%
  tt()
pairs Df SumsOfSqs F.Model R2 p.value p.adjusted sig
Control.1_Acclimation vs Treatment.1_Acclimation 1 0.2759858 0.9928976 0.06208366 0.488 1.000
Control.1_Acclimation vs Hot_control.1_Acclimation 1 1.4146464 6.5078325 0.30257965 0.001 0.015 .
Control.1_Acclimation vs Treatment.6_Post-FMT2 1 1.1524353 4.9536068 0.24825621 0.001 0.015 .
Control.1_Acclimation vs Control.6_Post-FMT2 1 0.4748999 1.9609749 0.11561688 0.005 0.075
Control.1_Acclimation vs Hot_control.6_Post-FMT2 1 1.3503168 5.4081420 0.26499923 0.001 0.015 .
Treatment.1_Acclimation vs Hot_control.1_Acclimation 1 1.6347704 6.8326887 0.29925029 0.001 0.015 .
Treatment.1_Acclimation vs Treatment.6_Post-FMT2 1 1.3540292 5.3398081 0.25022756 0.001 0.015 .
Treatment.1_Acclimation vs Control.6_Post-FMT2 1 0.6311089 2.4041625 0.13063146 0.001 0.015 .
Treatment.1_Acclimation vs Hot_control.6_Post-FMT2 1 1.6125755 5.9825981 0.27215155 0.001 0.015 .
Hot_control.1_Acclimation vs Treatment.6_Post-FMT2 1 0.6202327 3.1519868 0.16457754 0.001 0.015 .
Hot_control.1_Acclimation vs Control.6_Post-FMT2 1 1.5701179 7.6327037 0.32297209 0.001 0.015 .
Hot_control.1_Acclimation vs Hot_control.6_Post-FMT2 1 0.3634438 1.7083388 0.09647087 0.042 0.630
Treatment.6_Post-FMT2 vs Control.6_Post-FMT2 1 1.0227481 4.6483346 0.22511910 0.001 0.015 .
Treatment.6_Post-FMT2 vs Hot_control.6_Post-FMT2 1 0.5010202 2.2065321 0.12119453 0.002 0.030 .
Control.6_Post-FMT2 vs Hot_control.6_Post-FMT2 1 1.3619424 5.7710313 0.26507845 0.001 0.015 .
6.3.4.2.4 Phylogenetic
phylo_post4 <- as.matrix(beta_q1p$S)
phylo_post4 <- as.dist(phylo_post4[rownames(phylo_post4) %in% samples_to_keep_post4,
               colnames(phylo_post4) %in% samples_to_keep_post4])
betadisper(phylo_post4, subset_meta_post4$Population) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq  Mean Sq     F N.Perm Pr(>F)  
Groups     1 0.03098 0.030984 2.885    999   0.08 .
Residuals 51 0.54772 0.010740                      
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
         Cold_wet Hot_dry
Cold_wet            0.088
Hot_dry    0.0955        
adonis2(phylo_post4 ~ Population*time_point,
        data = subset_meta_post4 %>% arrange(match(Tube_code,labels(phylo_post4))),
        permutations = 999,
        strata = subset_meta_post4 %>% arrange(match(Tube_code,labels(phylo_post4))) %>% pull(individual)) %>%
        broom::tidy() %>%
        tt()
term df SumOfSqs R2 statistic p.value
Model 3 0.6466392 0.2442162 5.277784 0.001
Residual 49 2.0011759 0.7557838 NA NA
Total 52 2.6478151 1.0000000 NA NA
pairwise <- pairwise.adonis(phylo_post4, subset_meta_post4_arrange$type_time, perm=999)
pairwise%>%
  tt()
pairs Df SumsOfSqs F.Model R2 p.value p.adjusted sig
Control.1_Acclimation vs Treatment.1_Acclimation 1 0.05690611 0.789330 0.04999136 0.510 1.000
Control.1_Acclimation vs Hot_control.1_Acclimation 1 0.10866997 3.054731 0.16919285 0.012 0.180
Control.1_Acclimation vs Treatment.6_Post-FMT2 1 0.15591805 4.379209 0.22597458 0.007 0.105
Control.1_Acclimation vs Control.6_Post-FMT2 1 0.07099742 1.879364 0.11134092 0.123 1.000
Control.1_Acclimation vs Hot_control.6_Post-FMT2 1 0.18682367 4.878754 0.24542555 0.001 0.015 .
Treatment.1_Acclimation vs Hot_control.1_Acclimation 1 0.23108846 4.052184 0.20208192 0.003 0.045 .
Treatment.1_Acclimation vs Treatment.6_Post-FMT2 1 0.36496892 6.396667 0.28560797 0.001 0.015 .
Treatment.1_Acclimation vs Control.6_Post-FMT2 1 0.22628210 3.829222 0.19311005 0.017 0.255
Treatment.1_Acclimation vs Hot_control.6_Post-FMT2 1 0.34830814 5.846334 0.26761166 0.001 0.015 .
Hot_control.1_Acclimation vs Treatment.6_Post-FMT2 1 0.10002871 4.383624 0.21505615 0.001 0.015 .
Hot_control.1_Acclimation vs Control.6_Post-FMT2 1 0.12577510 5.060129 0.24027055 0.001 0.015 .
Hot_control.1_Acclimation vs Hot_control.6_Post-FMT2 1 0.06334378 2.499774 0.13512455 0.024 0.360
Treatment.6_Post-FMT2 vs Control.6_Post-FMT2 1 0.05927454 2.382025 0.12958449 0.024 0.360
Treatment.6_Post-FMT2 vs Hot_control.6_Post-FMT2 1 0.06906280 2.722460 0.14541146 0.006 0.090
Control.6_Post-FMT2 vs Hot_control.6_Post-FMT2 1 0.11081709 4.043656 0.20174244 0.002 0.030 .
6.3.4.2.5 Functional
func_post4 <- as.matrix(beta_q1f$S)
func_post4 <- as.dist(func_post4[rownames(func_post4) %in% samples_to_keep_post4,
               colnames(func_post4) %in% samples_to_keep_post4])
betadisper(func_post4, subset_meta_post4$Population) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq   Mean Sq     F N.Perm Pr(>F)
Groups     1 0.00002 0.0000207 0.001    999  0.978
Residuals 51 1.04916 0.0205717                    

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
         Cold_wet Hot_dry
Cold_wet            0.979
Hot_dry    0.9748        
adonis2(func_post4 ~ Population*time_point,
        data = subset_meta_post4 %>% arrange(match(Tube_code,labels(func_post4))),
        permutations = 999,
        strata = subset_meta_post4 %>% arrange(match(Tube_code,labels(func_post4))) %>% pull(individual)) %>%
        broom::tidy() %>%
        tt()
term df SumOfSqs R2 statistic p.value
Model 3 0.02728193 0.0197108 0.3284164 0.569
Residual 49 1.35682891 0.9802892 NA NA
Total 52 1.38411084 1.0000000 NA NA
pairwise <- pairwise.adonis(func_post4, subset_meta_post4_arrange$type_time, perm=999)
pairwise%>%
  tt()
pairs Df SumsOfSqs F.Model R2 p.value p.adjusted sig
Control.1_Acclimation vs Treatment.1_Acclimation 1 2.213850e-02 0.64325180 0.0411200822 0.443 1
Control.1_Acclimation vs Hot_control.1_Acclimation 1 1.746742e-02 0.54472723 0.0350425723 0.427 1
Control.1_Acclimation vs Treatment.6_Post-FMT2 1 3.674404e-02 1.08497265 0.0674525642 0.322 1
Control.1_Acclimation vs Control.6_Post-FMT2 1 2.983798e-02 0.93187248 0.0584910835 0.354 1
Control.1_Acclimation vs Hot_control.6_Post-FMT2 1 1.407466e-02 0.28346685 0.0185472873 0.517 1
Treatment.1_Acclimation vs Hot_control.1_Acclimation 1 1.373886e-03 0.07238159 0.0045034760 0.645 1
Treatment.1_Acclimation vs Treatment.6_Post-FMT2 1 -4.929805e-04 -0.02385162 -0.0014929516 0.713 1
Treatment.1_Acclimation vs Control.6_Post-FMT2 1 2.443617e-03 0.12903840 0.0080003779 0.576 1
Treatment.1_Acclimation vs Hot_control.6_Post-FMT2 1 6.968380e-03 0.19647180 0.0121305310 0.576 1
Hot_control.1_Acclimation vs Treatment.6_Post-FMT2 1 3.866607e-04 0.02093980 0.0013070267 0.673 1
Hot_control.1_Acclimation vs Control.6_Post-FMT2 1 -5.633463e-05 -0.00336651 -0.0002104511 0.681 1
Hot_control.1_Acclimation vs Hot_control.6_Post-FMT2 1 2.011448e-03 0.06046867 0.0037650628 0.734 1
Treatment.6_Post-FMT2 vs Control.6_Post-FMT2 1 -8.527330e-03 -0.46290555 -0.0297935723 0.849 1
Treatment.6_Post-FMT2 vs Hot_control.6_Post-FMT2 1 -1.648721e-03 -0.04717131 -0.0029569243 0.890 1
Control.6_Post-FMT2 vs Hot_control.6_Post-FMT2 1 4.367477e-03 0.13147026 0.0081499244 0.707 1
beta_richness_nmds_post4 <- richness_post4 %>%
  metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
  scores() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(subset_meta_post4, by = c("sample" = "Tube_code"))

beta_neutral_nmds_post4 <- neutral_post4 %>%
  metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
  scores() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(subset_meta_post4, by = c("sample" = "Tube_code"))

beta_phylo_nmds_post4 <- phylo_post4 %>%
  metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
  scores() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(subset_meta_post4, by = join_by(sample == Tube_code))

beta_func_nmds_post4 <- func_post4 %>%
  metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
  scores() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(subset_meta_post4, by = join_by(sample == Tube_code))

6.3.4.3 Comparison between Acclimation vs Post-FMT1 and Post-FMT2

6.3.4.3.1 Number of samples used
[1] 79
6.3.4.3.2 Richness
richness_post6 <- as.matrix(beta_q0n$S)
richness_post6 <- as.dist(richness_post6[rownames(richness_post6) %in% samples_to_keep_post6,
               colnames(richness_post6) %in% samples_to_keep_post6])
betadisper(richness_post6, subset_meta_post6$Population) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq  Mean Sq      F N.Perm Pr(>F)    
Groups     1 0.09792 0.097923 18.607    999  0.001 ***
Residuals 77 0.40523 0.005263                         
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
          Cold_wet Hot_dry
Cold_wet             0.001
Hot_dry  4.714e-05        
adonis2(richness_post6 ~ type*time_point,
        data = subset_meta_post6 %>% arrange(match(Tube_code,labels(richness_post6))),
        permutations = 999,
        strata = subset_meta_post6 %>% arrange(match(Tube_code,labels(richness_post6))) %>% pull(individual)) %>%
        broom::tidy() %>%
        tt()
term df SumOfSqs R2 statistic p.value
Model 8 6.395467 0.2451322 2.841433 0.001
Residual 70 19.694403 0.7548678 NA NA
Total 78 26.089870 1.0000000 NA NA
#Arrange of metadata dataframe
subset_meta_post6_arrange <- column_to_rownames(subset_meta_post6, "Tube_code")
subset_meta_post6_arrange<-subset_meta_post6_arrange[labels(richness_post6),]

pairwise <- pairwise.adonis(richness_post6, subset_meta_post6_arrange$type_time, perm=999)
pairwise%>%
  tt()
pairs Df SumsOfSqs F.Model R2 p.value p.adjusted sig
Control.1_Acclimation vs Treatment.1_Acclimation 1 0.3657243 1.123239 0.06966584 0.218 1.000
Control.1_Acclimation vs Hot_control.1_Acclimation 1 1.2605860 4.943717 0.24788342 0.001 0.036 .
Control.1_Acclimation vs Control.5_Post-FMT1 1 0.6804988 2.100611 0.12283837 0.002 0.072
Control.1_Acclimation vs Treatment.5_Post-FMT1 1 0.8481258 2.673290 0.16033371 0.001 0.036 .
Control.1_Acclimation vs Hot_control.5_Post-FMT1 1 1.1180661 3.809275 0.20252111 0.001 0.036 .
Control.1_Acclimation vs Treatment.6_Post-FMT2 1 0.8072604 2.940901 0.16392161 0.001 0.036 .
Control.1_Acclimation vs Control.6_Post-FMT2 1 0.4817387 1.660775 0.09968176 0.027 0.972
Control.1_Acclimation vs Hot_control.6_Post-FMT2 1 1.1179704 3.885459 0.20573812 0.001 0.036 .
Treatment.1_Acclimation vs Hot_control.1_Acclimation 1 1.3606630 5.087152 0.24124415 0.002 0.072
Treatment.1_Acclimation vs Control.5_Post-FMT1 1 0.7216200 2.172734 0.11956009 0.001 0.036 .
Treatment.1_Acclimation vs Treatment.5_Post-FMT1 1 0.9551308 2.926054 0.16322910 0.001 0.036 .
Treatment.1_Acclimation vs Hot_control.5_Post-FMT1 1 1.2263345 4.039487 0.20157637 0.001 0.036 .
Treatment.1_Acclimation vs Treatment.6_Post-FMT2 1 0.9130048 3.195028 0.16645080 0.001 0.036 .
Treatment.1_Acclimation vs Control.6_Post-FMT2 1 0.5959230 1.984036 0.11032208 0.002 0.072
Treatment.1_Acclimation vs Hot_control.6_Post-FMT2 1 1.2747787 4.275366 0.21086503 0.001 0.036 .
Hot_control.1_Acclimation vs Control.5_Post-FMT1 1 1.4319792 5.384836 0.25180628 0.001 0.036 .
Hot_control.1_Acclimation vs Treatment.5_Post-FMT1 1 0.8172413 3.194690 0.17558364 0.001 0.036 .
Hot_control.1_Acclimation vs Hot_control.5_Post-FMT1 1 0.5796135 2.441615 0.13239702 0.002 0.072
Hot_control.1_Acclimation vs Treatment.6_Post-FMT2 1 0.6397330 2.913695 0.15405213 0.001 0.036 .
Hot_control.1_Acclimation vs Control.6_Post-FMT2 1 1.4575447 6.224524 0.28007456 0.001 0.036 .
Hot_control.1_Acclimation vs Hot_control.6_Post-FMT2 1 0.3276169 1.412318 0.08111028 0.037 1.000
Control.5_Post-FMT1 vs Treatment.5_Post-FMT1 1 0.5615418 1.729004 0.10335366 0.017 0.612
Control.5_Post-FMT1 vs Hot_control.5_Post-FMT1 1 0.8438429 2.793772 0.14865413 0.002 0.072
Control.5_Post-FMT1 vs Treatment.6_Post-FMT2 1 0.7628135 2.683925 0.14364890 0.001 0.036 .
Control.5_Post-FMT1 vs Control.6_Post-FMT2 1 0.3432605 1.148733 0.06698647 0.246 1.000
Control.5_Post-FMT1 vs Hot_control.6_Post-FMT2 1 1.1269580 3.799256 0.19188884 0.001 0.036 .
Treatment.5_Post-FMT1 vs Hot_control.5_Post-FMT1 1 0.3734921 1.268929 0.07799710 0.115 1.000
Treatment.5_Post-FMT1 vs Treatment.6_Post-FMT2 1 0.3571397 1.297184 0.07959561 0.125 1.000
Treatment.5_Post-FMT1 vs Control.6_Post-FMT2 1 0.7769467 2.670898 0.15114670 0.001 0.036 .
Treatment.5_Post-FMT1 vs Hot_control.6_Post-FMT2 1 0.6502360 2.253407 0.13060650 0.001 0.036 .
Hot_control.5_Post-FMT1 vs Treatment.6_Post-FMT2 1 0.4132091 1.616138 0.09174188 0.009 0.324
Hot_control.5_Post-FMT1 vs Control.6_Post-FMT2 1 1.0163992 3.760571 0.19030682 0.001 0.036 .
Hot_control.5_Post-FMT1 vs Hot_control.6_Post-FMT2 1 0.2732563 1.019281 0.05988979 0.429 1.000
Treatment.6_Post-FMT2 vs Control.6_Post-FMT2 1 0.6463814 2.560441 0.13795154 0.002 0.072
Treatment.6_Post-FMT2 vs Hot_control.6_Post-FMT2 1 0.4796256 1.916520 0.10696943 0.004 0.144
Control.6_Post-FMT2 vs Hot_control.6_Post-FMT2 1 1.1305044 4.268317 0.21059061 0.001 0.036 .
6.3.4.3.3 Neutral
neutral_post6 <- as.matrix(beta_q1n$S)
neutral_post6 <- as.dist(neutral_post6[rownames(neutral_post6) %in% samples_to_keep_post6,
               colnames(neutral_post6) %in% samples_to_keep_post6])
betadisper(neutral_post6, subset_meta_post6$Population) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq  Mean Sq      F N.Perm Pr(>F)    
Groups     1 0.09524 0.095237 15.396    999  0.001 ***
Residuals 77 0.47632 0.006186                         
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
           Cold_wet Hot_dry
Cold_wet              0.001
Hot_dry  0.00018818        
adonis2(neutral_post6 ~ type*time_point,
        data = subset_meta_post6 %>% arrange(match(Tube_code,labels(neutral_post6))),
        permutations = 999,
        strata = subset_meta_post6 %>% arrange(match(Tube_code,labels(neutral_post6))) %>% pull(individual)) %>%
        broom::tidy() %>%
        tt()
term df SumOfSqs R2 statistic p.value
Model 8 7.736408 0.3121974 3.971673 0.001
Residual 70 17.044094 0.6878026 NA NA
Total 78 24.780501 1.0000000 NA NA
pairwise <- pairwise.adonis(neutral_post6, subset_meta_post6_arrange$type_time, perm=999)
pairwise%>%
  tt()
pairs Df SumsOfSqs F.Model R2 p.value p.adjusted sig
Control.1_Acclimation vs Treatment.1_Acclimation 1 0.2759858 0.9928976 0.06208366 0.482 1.000
Control.1_Acclimation vs Hot_control.1_Acclimation 1 1.4146464 6.5078325 0.30257965 0.001 0.036 .
Control.1_Acclimation vs Control.5_Post-FMT1 1 0.8153894 3.0970603 0.17113610 0.001 0.036 .
Control.1_Acclimation vs Treatment.5_Post-FMT1 1 1.1809241 4.4856470 0.24265567 0.002 0.072
Control.1_Acclimation vs Hot_control.5_Post-FMT1 1 1.4321524 5.7774260 0.27806264 0.001 0.036 .
Control.1_Acclimation vs Treatment.6_Post-FMT2 1 1.1524353 4.9536068 0.24825621 0.001 0.036 .
Control.1_Acclimation vs Control.6_Post-FMT2 1 0.4748999 1.9609749 0.11561688 0.008 0.288
Control.1_Acclimation vs Hot_control.6_Post-FMT2 1 1.3503168 5.4081420 0.26499923 0.001 0.036 .
Treatment.1_Acclimation vs Hot_control.1_Acclimation 1 1.6347704 6.8326887 0.29925029 0.001 0.036 .
Treatment.1_Acclimation vs Control.5_Post-FMT1 1 0.9517634 3.3715700 0.17404733 0.002 0.072
Treatment.1_Acclimation vs Treatment.5_Post-FMT1 1 1.3127773 4.6298256 0.23585668 0.001 0.036 .
Treatment.1_Acclimation vs Hot_control.5_Post-FMT1 1 1.6713369 6.2395460 0.28056085 0.001 0.036 .
Treatment.1_Acclimation vs Treatment.6_Post-FMT2 1 1.3540292 5.3398081 0.25022756 0.001 0.036 .
Treatment.1_Acclimation vs Control.6_Post-FMT2 1 0.6311089 2.4041625 0.13063146 0.001 0.036 .
Treatment.1_Acclimation vs Hot_control.6_Post-FMT2 1 1.6125755 5.9825981 0.27215155 0.001 0.036 .
Hot_control.1_Acclimation vs Control.5_Post-FMT1 1 1.5409781 6.8338056 0.29928456 0.001 0.036 .
Hot_control.1_Acclimation vs Treatment.5_Post-FMT1 1 0.9133614 4.0964534 0.21451383 0.001 0.036 .
Hot_control.1_Acclimation vs Hot_control.5_Post-FMT1 1 0.6954835 3.2951234 0.17077493 0.001 0.036 .
Hot_control.1_Acclimation vs Treatment.6_Post-FMT2 1 0.6202327 3.1519868 0.16457754 0.001 0.036 .
Hot_control.1_Acclimation vs Control.6_Post-FMT2 1 1.5701179 7.6327037 0.32297209 0.001 0.036 .
Hot_control.1_Acclimation vs Hot_control.6_Post-FMT2 1 0.3634438 1.7083388 0.09647087 0.029 1.000
Control.5_Post-FMT1 vs Treatment.5_Post-FMT1 1 0.6051778 2.2508491 0.13047758 0.016 0.576
Control.5_Post-FMT1 vs Hot_control.5_Post-FMT1 1 1.0528902 4.1436369 0.20570451 0.001 0.036 .
Control.5_Post-FMT1 vs Treatment.6_Post-FMT2 1 0.8908158 3.7146920 0.18842252 0.003 0.108
Control.5_Post-FMT1 vs Control.6_Post-FMT2 1 0.3860927 1.5521758 0.08843210 0.081 1.000
Control.5_Post-FMT1 vs Hot_control.6_Post-FMT2 1 1.3122237 5.1302726 0.24279254 0.001 0.036 .
Treatment.5_Post-FMT1 vs Hot_control.5_Post-FMT1 1 0.4150076 1.6372683 0.09840968 0.047 1.000
Treatment.5_Post-FMT1 vs Treatment.6_Post-FMT2 1 0.3157079 1.3252026 0.08117526 0.169 1.000
Treatment.5_Post-FMT1 vs Control.6_Post-FMT2 1 1.0579520 4.2700097 0.22158835 0.001 0.036 .
Treatment.5_Post-FMT1 vs Hot_control.6_Post-FMT2 1 0.7454015 2.9200493 0.16294873 0.002 0.072
Hot_control.5_Post-FMT1 vs Treatment.6_Post-FMT2 1 0.4377161 1.9421261 0.10824392 0.010 0.360
Hot_control.5_Post-FMT1 vs Control.6_Post-FMT2 1 1.3766597 5.8752789 0.26858075 0.001 0.036 .
Hot_control.5_Post-FMT1 vs Hot_control.6_Post-FMT2 1 0.3176516 1.3161367 0.07600637 0.147 1.000
Treatment.6_Post-FMT2 vs Control.6_Post-FMT2 1 1.0227481 4.6483346 0.22511910 0.001 0.036 .
Treatment.6_Post-FMT2 vs Hot_control.6_Post-FMT2 1 0.5010202 2.2065321 0.12119453 0.003 0.108
Control.6_Post-FMT2 vs Hot_control.6_Post-FMT2 1 1.3619424 5.7710313 0.26507845 0.001 0.036 .
6.3.4.3.4 Phylogenetic
phylo_post6 <- as.matrix(beta_q1p$S)
phylo_post6 <- as.dist(phylo_post6[rownames(phylo_post6) %in% samples_to_keep_post6,
               colnames(phylo_post6) %in% samples_to_keep_post6])
betadisper(phylo_post6, subset_meta_post6$Population) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq  Mean Sq      F N.Perm Pr(>F)
Groups     1 0.02554 0.025541 2.1633    999  0.165
Residuals 77 0.90909 0.011806                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
         Cold_wet Hot_dry
Cold_wet            0.165
Hot_dry   0.14542        
adonis2(phylo_post6 ~ type*time_point,
        data = subset_meta_post6 %>% arrange(match(Tube_code,labels(phylo_post6))),
        permutations = 999,
        strata = subset_meta_post6 %>% arrange(match(Tube_code,labels(phylo_post6))) %>% pull(individual)) %>%
        broom::tidy() %>%
        tt()
term df SumOfSqs R2 statistic p.value
Model 8 1.055866 0.2692287 3.223651 0.001
Residual 70 2.865952 0.7307713 NA NA
Total 78 3.921818 1.0000000 NA NA
pairwise <- pairwise.adonis(phylo_post6, subset_meta_post6_arrange$type_time, perm=999)
pairwise%>%
  tt()
pairs Df SumsOfSqs F.Model R2 p.value p.adjusted sig
Control.1_Acclimation vs Treatment.1_Acclimation 1 0.05690611 0.7893300 0.04999136 0.548 1.000
Control.1_Acclimation vs Hot_control.1_Acclimation 1 0.10866997 3.0547314 0.16919285 0.009 0.324
Control.1_Acclimation vs Control.5_Post-FMT1 1 0.11760287 2.9988032 0.16661126 0.024 0.864
Control.1_Acclimation vs Treatment.5_Post-FMT1 1 0.09624422 1.7247115 0.10968160 0.165 1.000
Control.1_Acclimation vs Hot_control.5_Post-FMT1 1 0.18586788 4.1904227 0.21836011 0.005 0.180
Control.1_Acclimation vs Treatment.6_Post-FMT2 1 0.15591805 4.3792085 0.22597458 0.007 0.252
Control.1_Acclimation vs Control.6_Post-FMT2 1 0.07099742 1.8793639 0.11134092 0.122 1.000
Control.1_Acclimation vs Hot_control.6_Post-FMT2 1 0.18682367 4.8787543 0.24542555 0.001 0.036 .
Treatment.1_Acclimation vs Hot_control.1_Acclimation 1 0.23108846 4.0521838 0.20208192 0.002 0.072
Treatment.1_Acclimation vs Control.5_Post-FMT1 1 0.26358465 4.3608960 0.21417997 0.009 0.324
Treatment.1_Acclimation vs Treatment.5_Post-FMT1 1 0.25319427 3.2738422 0.17915456 0.048 1.000
Treatment.1_Acclimation vs Hot_control.5_Post-FMT1 1 0.39050120 5.9837393 0.27218933 0.001 0.036 .
Treatment.1_Acclimation vs Treatment.6_Post-FMT2 1 0.36496892 6.3966666 0.28560797 0.001 0.036 .
Treatment.1_Acclimation vs Control.6_Post-FMT2 1 0.22628210 3.8292220 0.19311005 0.020 0.720
Treatment.1_Acclimation vs Hot_control.6_Post-FMT2 1 0.34830814 5.8463335 0.26761166 0.001 0.036 .
Hot_control.1_Acclimation vs Control.5_Post-FMT1 1 0.14203376 5.4200212 0.25303529 0.001 0.036 .
Hot_control.1_Acclimation vs Treatment.5_Post-FMT1 1 0.09666753 2.3682173 0.13635351 0.019 0.684
Hot_control.1_Acclimation vs Hot_control.5_Post-FMT1 1 0.09252600 2.9824958 0.15711821 0.005 0.180
Hot_control.1_Acclimation vs Treatment.6_Post-FMT2 1 0.10002871 4.3836237 0.21505615 0.001 0.036 .
Hot_control.1_Acclimation vs Control.6_Post-FMT2 1 0.12577510 5.0601287 0.24027055 0.001 0.036 .
Hot_control.1_Acclimation vs Hot_control.6_Post-FMT2 1 0.06334378 2.4997737 0.13512455 0.023 0.828
Control.5_Post-FMT1 vs Treatment.5_Post-FMT1 1 0.01842535 0.4144162 0.02688498 0.778 1.000
Control.5_Post-FMT1 vs Hot_control.5_Post-FMT1 1 0.05987967 1.7387847 0.09802164 0.120 1.000
Control.5_Post-FMT1 vs Treatment.6_Post-FMT2 1 0.07917244 3.0180046 0.15869197 0.009 0.324
Control.5_Post-FMT1 vs Control.6_Post-FMT2 1 0.04335491 1.5335604 0.08746429 0.189 1.000
Control.5_Post-FMT1 vs Hot_control.6_Post-FMT2 1 0.10783045 3.7500438 0.18987521 0.002 0.072
Treatment.5_Post-FMT1 vs Hot_control.5_Post-FMT1 1 0.03212966 0.6477782 0.04139746 0.682 1.000
Treatment.5_Post-FMT1 vs Treatment.6_Post-FMT2 1 0.06393539 1.5651817 0.09448624 0.138 1.000
Treatment.5_Post-FMT1 vs Control.6_Post-FMT2 1 0.05265949 1.2240203 0.07544494 0.287 1.000
Treatment.5_Post-FMT1 vs Hot_control.6_Post-FMT2 1 0.09753501 2.2402429 0.12994265 0.015 0.540
Hot_control.5_Post-FMT1 vs Treatment.6_Post-FMT2 1 0.07228545 2.3279593 0.12701683 0.037 1.000
Hot_control.5_Post-FMT1 vs Control.6_Post-FMT2 1 0.11759094 3.5538444 0.18174658 0.001 0.036 .
Hot_control.5_Post-FMT1 vs Hot_control.6_Post-FMT2 1 0.06667255 1.9859527 0.11041687 0.096 1.000
Treatment.6_Post-FMT2 vs Control.6_Post-FMT2 1 0.05927454 2.3820253 0.12958449 0.024 0.864
Treatment.6_Post-FMT2 vs Hot_control.6_Post-FMT2 1 0.06906280 2.7224602 0.14541146 0.004 0.144
Control.6_Post-FMT2 vs Hot_control.6_Post-FMT2 1 0.11081709 4.0436561 0.20174244 0.002 0.072
6.3.4.3.5 Functional
func_post6 <- as.matrix(beta_q1f$S)
func_post6 <- as.dist(func_post6[rownames(func_post6) %in% samples_to_keep_post6,
               colnames(func_post6) %in% samples_to_keep_post6])
betadisper(func_post6, subset_meta_post6$Population) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq   Mean Sq      F N.Perm Pr(>F)
Groups     1 0.00012 0.0001211 0.0064    999  0.931
Residuals 77 1.45678 0.0189193                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
         Cold_wet Hot_dry
Cold_wet            0.936
Hot_dry   0.93644        
adonis2(func_post6 ~ type*time_point,
        data = subset_meta_post6 %>% arrange(match(Tube_code,labels(func_post6))),
        permutations = 999,
        strata = subset_meta_post6 %>% arrange(match(Tube_code,labels(func_post6))) %>% pull(individual)) %>%
        broom::tidy() %>%
        tt()
term df SumOfSqs R2 statistic p.value
Model 8 0.1638096 0.07637627 0.7235547 0.563
Residual 70 1.9809613 0.92362373 NA NA
Total 78 2.1447709 1.00000000 NA NA
pairwise <- pairwise.adonis(func_post6, subset_meta_post6_arrange$type_time, perm=999)
pairwise%>%
  tt()
pairs Df SumsOfSqs F.Model R2 p.value p.adjusted sig
Control.1_Acclimation vs Treatment.1_Acclimation 1 2.213850e-02 0.6432517992 4.112008e-02 0.444 1
Control.1_Acclimation vs Hot_control.1_Acclimation 1 1.746742e-02 0.5447272282 3.504257e-02 0.409 1
Control.1_Acclimation vs Control.5_Post-FMT1 1 -7.536319e-03 -0.2066225492 -1.396723e-02 0.886 1
Control.1_Acclimation vs Treatment.5_Post-FMT1 1 1.283191e-01 3.4871170642 1.994106e-01 0.068 1
Control.1_Acclimation vs Hot_control.5_Post-FMT1 1 5.195054e-02 1.2654806328 7.780161e-02 0.291 1
Control.1_Acclimation vs Treatment.6_Post-FMT2 1 3.674404e-02 1.0849726511 6.745256e-02 0.315 1
Control.1_Acclimation vs Control.6_Post-FMT2 1 2.983798e-02 0.9318724842 5.849108e-02 0.385 1
Control.1_Acclimation vs Hot_control.6_Post-FMT2 1 1.407466e-02 0.2834668511 1.854729e-02 0.497 1
Treatment.1_Acclimation vs Hot_control.1_Acclimation 1 1.373886e-03 0.0723815851 4.503476e-03 0.621 1
Treatment.1_Acclimation vs Control.5_Post-FMT1 1 2.173211e-03 0.0940247450 5.842215e-03 0.600 1
Treatment.1_Acclimation vs Treatment.5_Post-FMT1 1 5.911966e-02 2.6246179279 1.489177e-01 0.162 1
Treatment.1_Acclimation vs Hot_control.5_Post-FMT1 1 1.091194e-02 0.3981698590 2.428136e-02 0.501 1
Treatment.1_Acclimation vs Treatment.6_Post-FMT2 1 -4.929805e-04 -0.0238516159 -1.492952e-03 0.714 1
Treatment.1_Acclimation vs Control.6_Post-FMT2 1 2.443617e-03 0.1290384018 8.000378e-03 0.592 1
Treatment.1_Acclimation vs Hot_control.6_Post-FMT2 1 6.968380e-03 0.1964718029 1.213053e-02 0.547 1
Hot_control.1_Acclimation vs Control.5_Post-FMT1 1 -2.303582e-03 -0.1101670946 -6.933181e-03 0.747 1
Hot_control.1_Acclimation vs Treatment.5_Post-FMT1 1 5.890814e-02 2.9198761676 1.629406e-01 0.147 1
Hot_control.1_Acclimation vs Hot_control.5_Post-FMT1 1 5.904278e-03 0.2342787606 1.443112e-02 0.548 1
Hot_control.1_Acclimation vs Treatment.6_Post-FMT2 1 3.866607e-04 0.0209397959 1.307027e-03 0.647 1
Hot_control.1_Acclimation vs Control.6_Post-FMT2 1 -5.633463e-05 -0.0033665098 -2.104511e-04 0.702 1
Hot_control.1_Acclimation vs Hot_control.6_Post-FMT2 1 2.011448e-03 0.0604686736 3.765063e-03 0.741 1
Control.5_Post-FMT1 vs Treatment.5_Post-FMT1 1 1.161466e-01 4.7247913225 2.395357e-01 0.078 1
Control.5_Post-FMT1 vs Hot_control.5_Post-FMT1 1 5.000930e-02 1.7048260738 9.629160e-02 0.240 1
Control.5_Post-FMT1 vs Treatment.6_Post-FMT2 1 2.292148e-02 1.0143427269 5.961692e-02 0.324 1
Control.5_Post-FMT1 vs Control.6_Post-FMT2 1 1.305120e-02 0.6254809273 3.762182e-02 0.417 1
Control.5_Post-FMT1 vs Hot_control.6_Post-FMT2 1 -8.269058e-03 -0.2211195871 -1.401364e-02 0.921 1
Treatment.5_Post-FMT1 vs Hot_control.5_Post-FMT1 1 1.235859e-02 0.4238120210 2.747777e-02 0.480 1
Treatment.5_Post-FMT1 vs Treatment.6_Post-FMT2 1 2.837514e-02 1.2912549771 7.926062e-02 0.288 1
Treatment.5_Post-FMT1 vs Control.6_Post-FMT2 1 3.882231e-02 1.9287804257 1.139350e-01 0.227 1
Treatment.5_Post-FMT1 vs Hot_control.6_Post-FMT2 1 5.671762e-02 1.5020399981 9.102147e-02 0.256 1
Hot_control.5_Post-FMT1 vs Treatment.6_Post-FMT2 1 -3.766327e-03 -0.1400672850 -8.831518e-03 0.725 1
Hot_control.5_Post-FMT1 vs Control.6_Post-FMT2 1 2.020599e-03 0.0803166527 4.994718e-03 0.618 1
Hot_control.5_Post-FMT1 vs Hot_control.6_Post-FMT2 1 5.871214e-06 0.0001408358 8.802161e-06 0.857 1
Treatment.6_Post-FMT2 vs Control.6_Post-FMT2 1 -8.527330e-03 -0.4629055475 -2.979357e-02 0.875 1
Treatment.6_Post-FMT2 vs Hot_control.6_Post-FMT2 1 -1.648721e-03 -0.0471713074 -2.956924e-03 0.900 1
Control.6_Post-FMT2 vs Hot_control.6_Post-FMT2 1 4.367477e-03 0.1314702630 8.149924e-03 0.710 1
beta_richness_nmds_post6 <- richness_post6 %>%
  metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
  scores() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(subset_meta_post6, by = c("sample" = "Tube_code"))

beta_neutral_nmds_post6 <- neutral_post6 %>%
  metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
  scores() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(subset_meta_post6, by = c("sample" = "Tube_code"))

beta_phylo_nmds_post6 <- phylo_post6 %>%
  metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
  scores() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(subset_meta_post6, by = join_by(sample == Tube_code))

beta_func_nmds_post6 <- func_post6 %>%
  metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
  scores() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(subset_meta_post6, by = join_by(sample == Tube_code))

6.3.5 5. Are there differences between the control and the treatment group?

6.3.5.1 After 1 week –> Post-FMT1

6.3.5.1.1 Number of samples used
[1] 26
6.3.5.1.2 Richness
richness_post1 <- as.matrix(beta_q0n$S)
richness_post1 <- as.dist(richness_post1[rownames(richness_post1) %in% samples_to_keep_post1,
               colnames(richness_post1) %in% samples_to_keep_post1])
betadisper(richness_post1, subset_meta_post1$type) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df   Sum Sq   Mean Sq      F N.Perm Pr(>F)
Groups     2 0.017675 0.0088373 2.3825    999  0.106
Residuals 23 0.085312 0.0037092                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
              Control Hot_control Treatment
Control                 0.0050000     0.652
Hot_control 0.0068795                 0.196
Treatment   0.6248469   0.2084296          
adonis2(richness_post1 ~ type,
        data = subset_meta_post1 %>% arrange(match(Tube_code,labels(richness_post1))),
        permutations = 999,
        strata = subset_meta_post1 %>% arrange(match(Tube_code,labels(richness_post1))) %>% pull(individual)) %>%
        broom::tidy() %>%
        tt()
term df SumOfSqs R2 statistic p.value
Model 2 1.195567 0.1448246 1.947534 1
Residual 23 7.059710 0.8551754 NA NA
Total 25 8.255277 1.0000000 NA NA
#Arrange of metadata dataframe
subset_meta_post1_arrange <- column_to_rownames(subset_meta_post1, "Tube_code")
subset_meta_post1_arrange<-subset_meta_post1_arrange[labels(richness_post1),]

pairwise <- pairwise.adonis(richness_post1, subset_meta_post1_arrange$type, perm=999)
pairwise%>%
  tt()
pairs Df SumsOfSqs F.Model R2 p.value p.adjusted sig
Control vs Treatment 1 0.5615418 1.729004 0.1033537 0.013 0.039 .
Control vs Hot_control 1 0.8438429 2.793772 0.1486541 0.001 0.003 *
Treatment vs Hot_control 1 0.3734921 1.268929 0.0779971 0.116 0.348
6.3.5.1.3 Neutral
neutral_post1 <- as.matrix(beta_q1n$S)
neutral_post1 <- as.dist(neutral_post1[rownames(neutral_post1) %in% samples_to_keep_post1,
               colnames(neutral_post1) %in% samples_to_keep_post1])
betadisper(neutral_post1, subset_meta_post1$type) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df   Sum Sq   Mean Sq      F N.Perm Pr(>F)
Groups     2 0.011001 0.0055005 0.6303    999  0.563
Residuals 23 0.200714 0.0087267                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
            Control Hot_control Treatment
Control                 0.22600     0.945
Hot_control 0.21166                 0.452
Treatment   0.95468     0.43604          
adonis2(neutral_post1 ~ type,
        data = subset_meta_post1 %>% arrange(match(Tube_code,labels(neutral_post1))),
        permutations = 999,
        strata = subset_meta_post1 %>% arrange(match(Tube_code,labels(neutral_post1))) %>% pull(individual)) %>%
        broom::tidy() %>%
        tt()
term df SumOfSqs R2 statistic p.value
Model 2 1.395968 0.1900228 2.697931 1
Residual 23 5.950350 0.8099772 NA NA
Total 25 7.346318 1.0000000 NA NA
pairwise <- pairwise.adonis(neutral_post1, subset_meta_post1_arrange$type, perm=999)
pairwise%>%
  tt()
pairs Df SumsOfSqs F.Model R2 p.value p.adjusted sig
Control vs Treatment 1 0.6051778 2.250849 0.13047758 0.012 0.036 .
Control vs Hot_control 1 1.0528902 4.143637 0.20570451 0.001 0.003 *
Treatment vs Hot_control 1 0.4150076 1.637268 0.09840968 0.053 0.159
6.3.5.1.4 Phylogenetic
phylo_post1 <- as.matrix(beta_q1p$S)
phylo_post1 <- as.dist(phylo_post1[rownames(phylo_post1) %in% samples_to_keep_post1,
               colnames(phylo_post1) %in% samples_to_keep_post1])
betadisper(phylo_post1, subset_meta_post1$type) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq   Mean Sq      F N.Perm Pr(>F)
Groups     2 0.00440 0.0021994 0.1369    999  0.898
Residuals 23 0.36941 0.0160614                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
            Control Hot_control Treatment
Control                 0.91500     0.690
Hot_control 0.91505                 0.788
Treatment   0.63312     0.73046          
adonis2(phylo_post1 ~ type,
        data = subset_meta_post1 %>% arrange(match(Tube_code,labels(phylo_post1))),
        permutations = 999,
        strata = subset_meta_post1 %>% arrange(match(Tube_code,labels(phylo_post1))) %>% pull(individual)) %>%
        broom::tidy() %>%
        tt()
term df SumOfSqs R2 statistic p.value
Model 2 0.07451036 0.07059466 0.8735033 1
Residual 23 0.98095695 0.92940534 NA NA
Total 25 1.05546731 1.00000000 NA NA
pairwise <- pairwise.adonis(phylo_post1, subset_meta_post1_arrange$type, perm=999)
pairwise%>%
  tt()
pairs Df SumsOfSqs F.Model R2 p.value p.adjusted sig
Control vs Treatment 1 0.01842535 0.4144162 0.02688498 0.792 1.00
Control vs Hot_control 1 0.05987967 1.7387847 0.09802164 0.130 0.39
Treatment vs Hot_control 1 0.03212966 0.6477782 0.04139746 0.701 1.00
6.3.5.1.5 Functional
func_post1 <- as.matrix(beta_q1f$S)
func_post1 <- as.dist(func_post1[rownames(func_post1) %in% samples_to_keep_post1,
               colnames(func_post1) %in% samples_to_keep_post1])
betadisper(func_post1, subset_meta_post1$type) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq   Mean Sq     F N.Perm Pr(>F)
Groups     2 0.00400 0.0020014 0.145    999  0.862
Residuals 23 0.31753 0.0138057                    

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
            Control Hot_control Treatment
Control                 0.59500     0.732
Hot_control 0.59817                 0.821
Treatment   0.75141     0.83718          
adonis2(func_post1 ~ type,
        data = subset_meta_post1 %>% arrange(match(Tube_code,labels(func_post1))),
        permutations = 999,
        strata = subset_meta_post1 %>% arrange(match(Tube_code,labels(func_post1))) %>% pull(individual)) %>%
        broom::tidy() %>%
        tt()
term df SumOfSqs R2 statistic p.value
Model 2 0.1186444 0.1568566 2.139435 1
Residual 23 0.6377435 0.8431434 NA NA
Total 25 0.7563879 1.0000000 NA NA
pairwise <- pairwise.adonis(func_post1, subset_meta_post1_arrange$type, perm=999)
pairwise%>%
  tt()
pairs Df SumsOfSqs F.Model R2 p.value p.adjusted sig
Control vs Treatment 1 0.11614656 4.724791 0.23953568 0.072 0.216
Control vs Hot_control 1 0.05000930 1.704826 0.09629160 0.228 0.684
Treatment vs Hot_control 1 0.01235859 0.423812 0.02747777 0.507 1.000
beta_richness_nmds_post1 <- richness_post1 %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(subset_meta_post1, by = join_by(sample == Tube_code))

beta_neutral_nmds_post1 <- neutral_post1 %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(subset_meta_post1, by = join_by(sample == Tube_code))

beta_phylogenetic_nmds_post1 <- phylo_post1 %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(subset_meta_post1, by = join_by(sample == Tube_code))

beta_functional_nmds_post1 <- func_post1 %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(subset_meta_post1, by = join_by(sample == Tube_code))
p0<-beta_richness_nmds_post1 %>%
            group_by(type) %>%
            mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
            mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
            ungroup() %>%
            ggplot(., aes(x=NMDS1,y=NMDS2, color=type)) +
        scale_color_manual(name="Type",
                       breaks=c("Control", "Hot_control", "Treatment"),
                       labels=c("Cold-Cold", "Hot-Hot", "Cold-Hot"),
                       values=c("#4477AA","#d57d2c","#76b183")) +
                geom_point(size=2) +
                geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
        labs(y = "NMDS2", x="NMDS1 \n Richness beta diversity") +
                theme_classic() +
                theme(legend.position="none")

p1<-beta_neutral_nmds_post1 %>%
            group_by(type) %>%
            mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
            mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
            ungroup() %>%
            ggplot(., aes(x=NMDS1,y=NMDS2, color=type)) +
        scale_color_manual(name="Type",
                       breaks=c("Control", "Hot_control", "Treatment"),
                       labels=c("Cold-Cold", "Hot-Hot", "Cold-Hot"),
                       values=c("#4477AA","#d57d2c","#76b183")) +
                geom_point(size=2) +
                geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
        labs(y = "NMDS2", x="NMDS1 \n Neutral beta diversity") +
                theme_classic() +
                theme(legend.position="none")
  
p2<-beta_phylogenetic_nmds_post1 %>%
            group_by(type) %>%
            mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
            mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
            ungroup() %>%
            ggplot(., aes(x=NMDS1,y=NMDS2, color=type)) +
                scale_color_manual(name="Type",
                       breaks=c("Control", "Hot_control", "Treatment"),
                       labels=c("Cold-Cold", "Hot-Hot", "Cold-Hot"),
                       values=c("#4477AA","#d57d2c","#76b183")) +
                geom_point(size=2) +
                geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
        labs(y= element_blank (), x="NMDS1 \n Phylogenetic beta diversity") +
                theme_classic() +
                theme(legend.position="none")

p3<-beta_functional_nmds_post1 %>%
            group_by(type) %>%
            mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
            mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
            ungroup() %>%
            ggplot(., aes(x=NMDS1,y=NMDS2, color=type)) +
                scale_color_manual(name="Type",
                       breaks=c("Control", "Hot_control", "Treatment"),
                       labels=c("Cold-Cold", "Hot-Hot", "Cold-Hot"),
                       values=c("#4477AA","#d57d2c","#76b183")) +
                geom_point(size=2) +
                geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
        labs(y= element_blank (), x="NMDS1 \n Functional beta diversity") +
                theme_classic()+
                theme(legend.position="none")
ggarrange(p0, p1, p2, p3, ncol=2, nrow=2, common.legend = TRUE, legend="right")

6.3.5.2 After 2 weeks –>Post-FMT2

post2 <- meta %>%
  filter(time_point == "6_Post-FMT2")

post2.counts <- temp_genome_counts[,which(colnames(temp_genome_counts) %in% rownames(post2))]
identical(sort(colnames(post2.counts)),sort(as.character(rownames(post2))))

post2_nmds <- sample_metadata %>%
  filter(time_point == "6_Post-FMT2")
6.3.5.2.1 Number of samples used
[1] 27
6.3.5.2.2 Richness
richness_post2 <- as.matrix(beta_q0n$S)
richness_post2 <- as.dist(richness_post2[rownames(richness_post2) %in% samples_to_keep_post2,
               colnames(richness_post2) %in% samples_to_keep_post2])
betadisper(richness_post2, subset_meta_post2$type) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df   Sum Sq   Mean Sq      F N.Perm Pr(>F)
Groups     2 0.002011 0.0010056 0.1982    999  0.845
Residuals 24 0.121775 0.0050740                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
            Control Hot_control Treatment
Control                 0.70000     0.827
Hot_control 0.67789                 0.633
Treatment   0.79246     0.59820          
adonis2(richness_post2 ~ type,
        data = subset_meta_post2 %>% arrange(match(Tube_code,labels(richness_post2))),
        permutations = 999,
        strata = subset_meta_post2 %>% arrange(match(Tube_code,labels(richness_post2))) %>% pull(individual)) %>%
        broom::tidy() %>%
        tt()
term df SumOfSqs R2 statistic p.value
Model 2 1.504341 0.1967776 2.939822 1
Residual 24 6.140538 0.8032224 NA NA
Total 26 7.644879 1.0000000 NA NA
#Arrange of metadata dataframe
subset_meta_post2_arrange <- column_to_rownames(subset_meta_post2, "Tube_code")
subset_meta_post2_arrange<-subset_meta_post2_arrange[labels(richness_post2),]

pairwise <- pairwise.adonis(richness_post2, subset_meta_post2_arrange$type, perm=999)
pairwise%>%
  tt()
pairs Df SumsOfSqs F.Model R2 p.value p.adjusted sig
Treatment vs Control 1 0.6463814 2.560441 0.1379515 0.001 0.003 *
Treatment vs Hot_control 1 0.4796256 1.916520 0.1069694 0.001 0.003 *
Control vs Hot_control 1 1.1305044 4.268317 0.2105906 0.001 0.003 *
6.3.5.2.3 Neutral
neutral_post2 <- as.matrix(beta_q1n$S)
neutral_post2 <- as.dist(neutral_post2[rownames(neutral_post2) %in% samples_to_keep_post2,
               colnames(neutral_post2) %in% samples_to_keep_post2])
betadisper(neutral_post2, subset_meta_post2$type) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df   Sum Sq   Mean Sq      F N.Perm Pr(>F)
Groups     2 0.008262 0.0041311 0.8024    999  0.458
Residuals 24 0.123559 0.0051483                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
            Control Hot_control Treatment
Control                 0.44000     0.649
Hot_control 0.44675                 0.249
Treatment   0.65989     0.25095          
adonis2(neutral_post2 ~ type,
        data = subset_meta_post2 %>% arrange(match(Tube_code,labels(neutral_post2))),
        permutations = 999,
        strata = subset_meta_post2 %>% arrange(match(Tube_code,labels(neutral_post2))) %>% pull(individual)) %>%
        broom::tidy() %>%
        tt()
term df SumOfSqs R2 statistic p.value
Model 2 1.923807 0.2603795 4.224537 1
Residual 24 5.464666 0.7396205 NA NA
Total 26 7.388473 1.0000000 NA NA
pairwise <- pairwise.adonis(neutral_post2, subset_meta_post2_arrange$type, perm=999)
pairwise%>%
  tt()
pairs Df SumsOfSqs F.Model R2 p.value p.adjusted sig
Treatment vs Control 1 1.0227481 4.648335 0.2251191 0.001 0.003 *
Treatment vs Hot_control 1 0.5010202 2.206532 0.1211945 0.003 0.009 *
Control vs Hot_control 1 1.3619424 5.771031 0.2650785 0.001 0.003 *
6.3.5.2.4 Phylogenetic
phylo_post2 <- as.matrix(beta_q1p$S)
phylo_post2 <- as.dist(phylo_post2[rownames(phylo_post2) %in% samples_to_keep_post2,
               colnames(phylo_post2) %in% samples_to_keep_post2])
betadisper(phylo_post2, subset_meta_post2$type) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df   Sum Sq   Mean Sq      F N.Perm Pr(>F)
Groups     2 0.000407 0.0002034 0.0487    999  0.947
Residuals 24 0.100305 0.0041794                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
            Control Hot_control Treatment
Control                 0.92800     0.846
Hot_control 0.93765                 0.737
Treatment   0.83933     0.76015          
adonis2(phylo_post2 ~ type,
        data = subset_meta_post2 %>% arrange(match(Tube_code,labels(phylo_post2))),
        permutations = 999,
        strata = subset_meta_post2 %>% arrange(match(Tube_code,labels(phylo_post2))) %>% pull(individual)) %>%
        broom::tidy() %>%
        tt()
term df SumOfSqs R2 statistic p.value
Model 2 0.1594363 0.2042241 3.079623 1
Residual 24 0.6212564 0.7957759 NA NA
Total 26 0.7806927 1.0000000 NA NA
pairwise <- pairwise.adonis(phylo_post2, subset_meta_post2_arrange$type, perm=999)
pairwise%>%
  tt()
pairs Df SumsOfSqs F.Model R2 p.value p.adjusted sig
Treatment vs Control 1 0.05927454 2.382025 0.1295845 0.027 0.081
Treatment vs Hot_control 1 0.06906280 2.722460 0.1454115 0.004 0.012 .
Control vs Hot_control 1 0.11081709 4.043656 0.2017424 0.002 0.006 *
6.3.5.2.5 Functional
func_post2 <- as.matrix(beta_q1f$S)
func_post2 <- as.dist(func_post2[rownames(func_post2) %in% samples_to_keep_post2,
               colnames(func_post2) %in% samples_to_keep_post2])
betadisper(func_post2, subset_meta_post2$type) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq   Mean Sq      F N.Perm Pr(>F)
Groups     2 0.01126 0.0056302 0.2861    999  0.813
Residuals 24 0.47233 0.0196806                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
            Control Hot_control Treatment
Control                 0.55500     0.681
Hot_control 0.48255                 0.816
Treatment   0.60116     0.75643          
adonis2(func_post2 ~ type,
        data = subset_meta_post2 %>% arrange(match(Tube_code,labels(func_post2))),
        permutations = 999,
        strata = subset_meta_post2 %>% arrange(match(Tube_code,labels(func_post2))) %>% pull(individual)) %>%
        broom::tidy() %>%
        tt()
term df SumOfSqs R2 statistic p.value
Model 2 -0.003872383 -0.005621319 -0.06707875 1
Residual 24 0.692746821 1.005621319 NA NA
Total 26 0.688874438 1.000000000 NA NA
pairwise <- pairwise.adonis(func_post2, subset_meta_post2_arrange$type, perm=999)
pairwise%>%
  tt()
pairs Df SumsOfSqs F.Model R2 p.value p.adjusted sig
Treatment vs Control 1 -0.008527330 -0.46290555 -0.029793572 0.836 1
Treatment vs Hot_control 1 -0.001648721 -0.04717131 -0.002956924 0.906 1
Control vs Hot_control 1 0.004367477 0.13147026 0.008149924 0.699 1
beta_richness_nmds_post2 <- richness_post2 %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(subset_meta_post2, by = join_by(sample == Tube_code))

beta_neutral_nmds_post2 <- neutral_post2 %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(subset_meta_post2, by = join_by(sample == Tube_code))

beta_phylogenetic_nmds_post2 <- phylo_post2 %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(subset_meta_post2, by = join_by(sample == Tube_code))

beta_functional_nmds_post2 <- func_post2 %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(subset_meta_post2, by = join_by(sample == Tube_code))
p0<-beta_richness_nmds_post2 %>%
            group_by(type) %>%
            mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
            mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
            ungroup() %>%
            ggplot(., aes(x=NMDS1,y=NMDS2, color=type)) +
        scale_color_manual(name="Type",
                       breaks=c("Control", "Hot_control", "Treatment"),
                       labels=c("Cold-Cold", "Hot-Hot", "Cold-Hot"),
                       values=c("#4477AA","#d57d2c","#76b183")) +
                geom_point(size=2) +
                geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
        labs(y = "NMDS2", x="NMDS1 \n Richness beta diversity") +
                theme_classic() +
                theme(legend.position="none")

p1<-beta_neutral_nmds_post2 %>%
            group_by(type) %>%
            mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
            mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
            ungroup() %>%
            ggplot(., aes(x=NMDS1,y=NMDS2, color=type)) +
        scale_color_manual(name="Type",
                       breaks=c("Control", "Hot_control", "Treatment"),
                       labels=c("Cold-Cold", "Hot-Hot", "Cold-Hot"),
                       values=c("#4477AA","#d57d2c","#76b183")) +
                geom_point(size=2) +
                geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
        labs(y = "NMDS2", x="NMDS1 \n Neutral beta diversity") +
                theme_classic() +
                theme(legend.position="none")
  
p2<-beta_phylogenetic_nmds_post2 %>%
            group_by(type) %>%
            mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
            mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
            ungroup() %>%
            ggplot(., aes(x=NMDS1,y=NMDS2, color=type)) +
                scale_color_manual(name="Type",
                       breaks=c("Control", "Hot_control", "Treatment"),
                       labels=c("Cold-Cold", "Hot-Hot", "Cold-Hot"),
                       values=c("#4477AA","#d57d2c","#76b183")) +
                geom_point(size=2) +
                geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
        labs(y= element_blank (), x="NMDS1 \n Phylogenetic beta diversity") +
                theme_classic() +
                theme(legend.position="none")

p3<-beta_functional_nmds_post2 %>%
            group_by(type) %>%
            mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
            mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
            ungroup() %>%
            ggplot(., aes(x=NMDS1,y=NMDS2, color=type)) +
                scale_color_manual(name="Type",
                       breaks=c("Control", "Hot_control", "Treatment"),
                       labels=c("Cold-Cold", "Hot-Hot", "Cold-Hot"),
                       values=c("#4477AA","#d57d2c","#76b183")) +
                geom_point(size=2) +
                geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
        labs(y= element_blank (), x="NMDS1 \n Functional beta diversity") +
                theme_classic()+
                theme(legend.position="none")
ggarrange(p0, p1, p2, p3, ncol=2, nrow=2, common.legend = TRUE, legend="right")

6.3.5.3 Post1 vs Post2

post5 <- meta %>%
  filter(time_point == "6_Post-FMT2" | time_point == "5_Post-FMT1")

post5.counts <- temp_genome_counts[,which(colnames(temp_genome_counts) %in% rownames(post5))]
identical(sort(colnames(post5.counts)),sort(as.character(rownames(post5))))

post5_nmds <- sample_metadata %>%
  filter(time_point == "6_Post-FMT2"| time_point == "5_Post-FMT1")
6.3.5.3.1 Number of samples used
[1] 53
6.3.5.3.2 Richness
richness_post5 <- as.matrix(beta_q0n$S)
richness_post5 <- as.dist(richness_post5[rownames(richness_post5) %in% samples_to_keep_post5,
               colnames(richness_post5) %in% samples_to_keep_post5])
betadisper(richness_post5, subset_meta_post5$type) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq   Mean Sq      F N.Perm Pr(>F)
Groups     2 0.01841 0.0092048 1.7364    999  0.178
Residuals 50 0.26505 0.0053010                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
             Control Hot_control Treatment
Control                 0.033000     0.714
Hot_control 0.039117                 0.229
Treatment   0.716358    0.218648          
adonis2(richness_post5 ~ type*time_point,
        data = subset_meta_post5 %>% arrange(match(Tube_code,labels(richness_post5))),
        permutations = 999,
        strata = subset_meta_post5 %>% arrange(match(Tube_code,labels(richness_post5))) %>% pull(individual)) %>%
        broom::tidy() %>%
        tt()
term df SumOfSqs R2 statistic p.value
Model 5 3.261916 0.1981463 2.322836 0.001
Residual 47 13.200248 0.8018537 NA NA
Total 52 16.462164 1.0000000 NA NA
#Arrange of metadata dataframe
subset_meta_post5_arrange <- column_to_rownames(subset_meta_post5, "Tube_code")
subset_meta_post5_arrange<-subset_meta_post5_arrange[labels(richness_post5),]

pairwise <- pairwise.adonis(richness_post5, subset_meta_post5_arrange$type_time, perm=999)
pairwise%>%
  tt()
pairs Df SumsOfSqs F.Model R2 p.value p.adjusted sig
Control.5_Post-FMT1 vs Treatment.5_Post-FMT1 1 0.5615418 1.729004 0.10335366 0.020 0.300
Control.5_Post-FMT1 vs Hot_control.5_Post-FMT1 1 0.8438429 2.793772 0.14865413 0.001 0.015 .
Control.5_Post-FMT1 vs Treatment.6_Post-FMT2 1 0.7628135 2.683925 0.14364890 0.001 0.015 .
Control.5_Post-FMT1 vs Control.6_Post-FMT2 1 0.3432605 1.148733 0.06698647 0.274 1.000
Control.5_Post-FMT1 vs Hot_control.6_Post-FMT2 1 1.1269580 3.799256 0.19188884 0.002 0.030 .
Treatment.5_Post-FMT1 vs Hot_control.5_Post-FMT1 1 0.3734921 1.268929 0.07799710 0.117 1.000
Treatment.5_Post-FMT1 vs Treatment.6_Post-FMT2 1 0.3571397 1.297184 0.07959561 0.128 1.000
Treatment.5_Post-FMT1 vs Control.6_Post-FMT2 1 0.7769467 2.670898 0.15114670 0.001 0.015 .
Treatment.5_Post-FMT1 vs Hot_control.6_Post-FMT2 1 0.6502360 2.253407 0.13060650 0.001 0.015 .
Hot_control.5_Post-FMT1 vs Treatment.6_Post-FMT2 1 0.4132091 1.616138 0.09174188 0.008 0.120
Hot_control.5_Post-FMT1 vs Control.6_Post-FMT2 1 1.0163992 3.760571 0.19030682 0.001 0.015 .
Hot_control.5_Post-FMT1 vs Hot_control.6_Post-FMT2 1 0.2732563 1.019281 0.05988979 0.413 1.000
Treatment.6_Post-FMT2 vs Control.6_Post-FMT2 1 0.6463814 2.560441 0.13795154 0.001 0.015 .
Treatment.6_Post-FMT2 vs Hot_control.6_Post-FMT2 1 0.4796256 1.916520 0.10696943 0.001 0.015 .
Control.6_Post-FMT2 vs Hot_control.6_Post-FMT2 1 1.1305044 4.268317 0.21059061 0.001 0.015 .
6.3.5.3.3 Neutral
neutral_post5 <- as.matrix(beta_q1n$S)
neutral_post5 <- as.dist(neutral_post5[rownames(neutral_post5) %in% samples_to_keep_post5,
               colnames(neutral_post5) %in% samples_to_keep_post5])
betadisper(neutral_post5, subset_meta_post5$type) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq   Mean Sq     F N.Perm Pr(>F)
Groups     2 0.01992 0.0099587 1.565    999  0.202
Residuals 50 0.31818 0.0063636                    

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
            Control Hot_control Treatment
Control                 0.10400     0.881
Hot_control 0.10701                 0.180
Treatment   0.87156     0.17449          
adonis2(neutral_post5 ~ type*time_point,
        data = subset_meta_post5 %>% arrange(match(Tube_code,labels(neutral_post5))),
        permutations = 999,
        strata = subset_meta_post5 %>% arrange(match(Tube_code,labels(neutral_post5))) %>% pull(individual)) %>%
        broom::tidy() %>%
        tt()
term df SumOfSqs R2 statistic p.value
Model 5 3.947979 0.2569798 3.251069 0.001
Residual 47 11.415016 0.7430202 NA NA
Total 52 15.362995 1.0000000 NA NA
pairwise <- pairwise.adonis(neutral_post5, subset_meta_post5_arrange$type_time, perm=999)
pairwise%>%
  tt()
pairs Df SumsOfSqs F.Model R2 p.value p.adjusted sig
Control.5_Post-FMT1 vs Treatment.5_Post-FMT1 1 0.6051778 2.250849 0.13047758 0.015 0.225
Control.5_Post-FMT1 vs Hot_control.5_Post-FMT1 1 1.0528902 4.143637 0.20570451 0.001 0.015 .
Control.5_Post-FMT1 vs Treatment.6_Post-FMT2 1 0.8908158 3.714692 0.18842252 0.002 0.030 .
Control.5_Post-FMT1 vs Control.6_Post-FMT2 1 0.3860927 1.552176 0.08843210 0.077 1.000
Control.5_Post-FMT1 vs Hot_control.6_Post-FMT2 1 1.3122237 5.130273 0.24279254 0.001 0.015 .
Treatment.5_Post-FMT1 vs Hot_control.5_Post-FMT1 1 0.4150076 1.637268 0.09840968 0.054 0.810
Treatment.5_Post-FMT1 vs Treatment.6_Post-FMT2 1 0.3157079 1.325203 0.08117526 0.158 1.000
Treatment.5_Post-FMT1 vs Control.6_Post-FMT2 1 1.0579520 4.270010 0.22158835 0.001 0.015 .
Treatment.5_Post-FMT1 vs Hot_control.6_Post-FMT2 1 0.7454015 2.920049 0.16294873 0.001 0.015 .
Hot_control.5_Post-FMT1 vs Treatment.6_Post-FMT2 1 0.4377161 1.942126 0.10824392 0.008 0.120
Hot_control.5_Post-FMT1 vs Control.6_Post-FMT2 1 1.3766597 5.875279 0.26858075 0.001 0.015 .
Hot_control.5_Post-FMT1 vs Hot_control.6_Post-FMT2 1 0.3176516 1.316137 0.07600637 0.182 1.000
Treatment.6_Post-FMT2 vs Control.6_Post-FMT2 1 1.0227481 4.648335 0.22511910 0.001 0.015 .
Treatment.6_Post-FMT2 vs Hot_control.6_Post-FMT2 1 0.5010202 2.206532 0.12119453 0.002 0.030 .
Control.6_Post-FMT2 vs Hot_control.6_Post-FMT2 1 1.3619424 5.771031 0.26507845 0.001 0.015 .
6.3.5.3.4 Phylogenetic
phylo_post5 <- as.matrix(beta_q1p$S)
phylo_post5 <- as.dist(phylo_post5[rownames(phylo_post5) %in% samples_to_keep_post5,
               colnames(phylo_post5) %in% samples_to_keep_post5])
betadisper(phylo_post5, subset_meta_post5$type) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq   Mean Sq      F N.Perm Pr(>F)
Groups     2 0.00051 0.0002543 0.0265    999   0.97
Residuals 50 0.47996 0.0095993                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
            Control Hot_control Treatment
Control                 0.87900     0.818
Hot_control 0.88926                 0.928
Treatment   0.82391     0.91902          
adonis2(phylo_post5 ~ type*time_point,
        data = subset_meta_post5 %>% arrange(match(Tube_code,labels(phylo_post5))),
        permutations = 999,
        strata = subset_meta_post5 %>% arrange(match(Tube_code,labels(phylo_post5))) %>% pull(individual)) %>%
        broom::tidy() %>%
        tt()
term df SumOfSqs R2 statistic p.value
Model 5 0.3518222 0.180049 2.0641 0.001
Residual 47 1.6022134 0.819951 NA NA
Total 52 1.9540356 1.000000 NA NA
pairwise <- pairwise.adonis(phylo_post5, subset_meta_post5_arrange$type_time, perm=999)
pairwise%>%
  tt()
pairs Df SumsOfSqs F.Model R2 p.value p.adjusted sig
Control.5_Post-FMT1 vs Treatment.5_Post-FMT1 1 0.01842535 0.4144162 0.02688498 0.806 1.000
Control.5_Post-FMT1 vs Hot_control.5_Post-FMT1 1 0.05987967 1.7387847 0.09802164 0.117 1.000
Control.5_Post-FMT1 vs Treatment.6_Post-FMT2 1 0.07917244 3.0180046 0.15869197 0.007 0.105
Control.5_Post-FMT1 vs Control.6_Post-FMT2 1 0.04335491 1.5335604 0.08746429 0.162 1.000
Control.5_Post-FMT1 vs Hot_control.6_Post-FMT2 1 0.10783045 3.7500438 0.18987521 0.001 0.015 .
Treatment.5_Post-FMT1 vs Hot_control.5_Post-FMT1 1 0.03212966 0.6477782 0.04139746 0.669 1.000
Treatment.5_Post-FMT1 vs Treatment.6_Post-FMT2 1 0.06393539 1.5651817 0.09448624 0.126 1.000
Treatment.5_Post-FMT1 vs Control.6_Post-FMT2 1 0.05265949 1.2240203 0.07544494 0.303 1.000
Treatment.5_Post-FMT1 vs Hot_control.6_Post-FMT2 1 0.09753501 2.2402429 0.12994265 0.010 0.150
Hot_control.5_Post-FMT1 vs Treatment.6_Post-FMT2 1 0.07228545 2.3279593 0.12701683 0.040 0.600
Hot_control.5_Post-FMT1 vs Control.6_Post-FMT2 1 0.11759094 3.5538444 0.18174658 0.001 0.015 .
Hot_control.5_Post-FMT1 vs Hot_control.6_Post-FMT2 1 0.06667255 1.9859527 0.11041687 0.090 1.000
Treatment.6_Post-FMT2 vs Control.6_Post-FMT2 1 0.05927454 2.3820253 0.12958449 0.022 0.330
Treatment.6_Post-FMT2 vs Hot_control.6_Post-FMT2 1 0.06906280 2.7224602 0.14541146 0.004 0.060
Control.6_Post-FMT2 vs Hot_control.6_Post-FMT2 1 0.11081709 4.0436561 0.20174244 0.002 0.030 .
6.3.5.3.5 Functional
func_post5 <- as.matrix(beta_q1f$S)
func_post5 <- as.dist(func_post5[rownames(func_post5) %in% samples_to_keep_post5,
               colnames(func_post5) %in% samples_to_keep_post5])
betadisper(func_post5, subset_meta_post5$type) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq   Mean Sq      F N.Perm Pr(>F)
Groups     2 0.00855 0.0042775 0.2522    999  0.787
Residuals 50 0.84815 0.0169631                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
            Control Hot_control Treatment
Control                 0.52500     0.611
Hot_control 0.50303                 0.844
Treatment   0.59083     0.82686          
adonis2(func_post5 ~ type*time_point,
        data = subset_meta_post5 %>% arrange(match(Tube_code,labels(func_post5))),
        permutations = 999,
        strata = subset_meta_post5 %>% arrange(match(Tube_code,labels(func_post5))) %>% pull(individual)) %>%
        broom::tidy() %>%
        tt()
term df SumOfSqs R2 statistic p.value
Model 5 0.1047944 0.07301299 0.7403794 0.471
Residual 47 1.3304903 0.92698701 NA NA
Total 52 1.4352848 1.00000000 NA NA
pairwise <- pairwise.adonis(func_post5, subset_meta_post5_arrange$type_time, perm=999)
pairwise%>%
  tt()
pairs Df SumsOfSqs F.Model R2 p.value p.adjusted sig
Control.5_Post-FMT1 vs Treatment.5_Post-FMT1 1 1.161466e-01 4.7247913225 2.395357e-01 0.075 1
Control.5_Post-FMT1 vs Hot_control.5_Post-FMT1 1 5.000930e-02 1.7048260738 9.629160e-02 0.230 1
Control.5_Post-FMT1 vs Treatment.6_Post-FMT2 1 2.292148e-02 1.0143427269 5.961692e-02 0.311 1
Control.5_Post-FMT1 vs Control.6_Post-FMT2 1 1.305120e-02 0.6254809273 3.762182e-02 0.404 1
Control.5_Post-FMT1 vs Hot_control.6_Post-FMT2 1 -8.269058e-03 -0.2211195871 -1.401364e-02 0.921 1
Treatment.5_Post-FMT1 vs Hot_control.5_Post-FMT1 1 1.235859e-02 0.4238120210 2.747777e-02 0.508 1
Treatment.5_Post-FMT1 vs Treatment.6_Post-FMT2 1 2.837514e-02 1.2912549771 7.926062e-02 0.301 1
Treatment.5_Post-FMT1 vs Control.6_Post-FMT2 1 3.882231e-02 1.9287804257 1.139350e-01 0.199 1
Treatment.5_Post-FMT1 vs Hot_control.6_Post-FMT2 1 5.671762e-02 1.5020399981 9.102147e-02 0.242 1
Hot_control.5_Post-FMT1 vs Treatment.6_Post-FMT2 1 -3.766327e-03 -0.1400672850 -8.831518e-03 0.731 1
Hot_control.5_Post-FMT1 vs Control.6_Post-FMT2 1 2.020599e-03 0.0803166527 4.994718e-03 0.607 1
Hot_control.5_Post-FMT1 vs Hot_control.6_Post-FMT2 1 5.871214e-06 0.0001408358 8.802161e-06 0.853 1
Treatment.6_Post-FMT2 vs Control.6_Post-FMT2 1 -8.527330e-03 -0.4629055475 -2.979357e-02 0.859 1
Treatment.6_Post-FMT2 vs Hot_control.6_Post-FMT2 1 -1.648721e-03 -0.0471713074 -2.956924e-03 0.910 1
Control.6_Post-FMT2 vs Hot_control.6_Post-FMT2 1 4.367477e-03 0.1314702630 8.149924e-03 0.696 1
beta_richness_nmds_post5 <- richness_post5 %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(subset_meta_post5, by = join_by(sample == Tube_code))

beta_neutral_nmds_post5 <- neutral_post5 %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(subset_meta_post5, by = join_by(sample == Tube_code))

beta_phylogenetic_nmds_post5 <- phylo_post5 %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(subset_meta_post5, by = join_by(sample == Tube_code))

beta_functional_nmds_post5 <- func_post5 %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(subset_meta_post5, by = join_by(sample == Tube_code))
p0<-beta_richness_nmds_post5 %>%
            group_by(type, time_point) %>%
            mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
            mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
            ungroup() %>%
            ggplot(., aes(x=NMDS1,y=NMDS2, color=type, shape=time_point)) +
        scale_color_manual(name="Type",
                       breaks=c("Control", "Hot_control", "Treatment"),
                       labels=c("Cold-Cold", "Hot-Hot", "Cold-Hot"),
                       values=c("#4477AA","#d57d2c","#76b183")) +
                geom_point(size=2) +
                geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
        labs(y = "NMDS2", x="NMDS1 \n Richness beta diversity") +
                theme_classic() +
                theme(legend.position="none")

p1<-beta_neutral_nmds_post5 %>%
            group_by(type, time_point) %>%
            mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
            mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
            ungroup() %>%
            ggplot(., aes(x=NMDS1,y=NMDS2, color=type, shape=time_point)) +
        scale_color_manual(name="Type",
                       breaks=c("Control", "Hot_control", "Treatment"),
                       labels=c("Cold-Cold", "Hot-Hot", "Cold-Hot"),
                       values=c("#4477AA","#d57d2c","#76b183")) +
                geom_point(size=2) +
                geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
        labs(y = "NMDS2", x="NMDS1 \n Neutral beta diversity") +
                theme_classic() +
                theme(legend.position="none")
  
p2<-beta_phylogenetic_nmds_post5 %>%
            group_by(type, time_point) %>%
            mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
            mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
            ungroup() %>%
            ggplot(., aes(x=NMDS1,y=NMDS2, color=type, shape=time_point)) +
                scale_color_manual(name="Type",
                       breaks=c("Control", "Hot_control", "Treatment"),
                       labels=c("Cold-Cold", "Hot-Hot", "Cold-Hot"),
                       values=c("#4477AA","#d57d2c","#76b183")) +
                geom_point(size=2) +
                geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
        labs(y= element_blank (), x="NMDS1 \n Phylogenetic beta diversity") +
                theme_classic() +
                theme(legend.position="none")

p3<-beta_functional_nmds_post5 %>%
            group_by(type, time_point) %>%
            mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
            mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
            ungroup() %>%
            ggplot(., aes(x=NMDS1,y=NMDS2, color=type, shape=time_point)) +
                scale_color_manual(name="Type",
                       breaks=c("Control", "Hot_control", "Treatment"),
                       labels=c("Cold-Cold", "Hot-Hot", "Cold-Hot"),
                       values=c("#4477AA","#d57d2c","#76b183")) +
                geom_point(size=2) +
                geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
        labs(y= element_blank (), x="NMDS1 \n Functional beta diversity") +
                theme_classic()+
                theme(legend.position="none")
ggarrange(p0, p1, p2, p3, ncol=2, nrow=2, common.legend = TRUE, legend="right")